Unified Lagrangian formulation for fluid and solid mechanics, fluid-structure interaction and coupled thermal problems using the PFEM

Abstract

toc El objectivo de la presente tesis es la derivación e implementación de una formulación unificada con elementos finitos para la solución de problemas de mecánica de fluidos y de sólidos, interacción fluido-estructura (Fluid-Structure Interaction (FSI)) y con acoplamiento térmico. El método unificado está basado en una formulación Lagrangiana estabilizada y las variables incçognitas son las velocidades y la presión. Cada paso de tiempo se soluciona a través de un esquema de dos pasos de tipo Gauss-Seidel. Primero se resuelven las ecuaciones de momento lineal por los incrementos de velocidad, luego se calculan las presiones en la configuración actualizada usando la ecuación de continuidad. Para los dominios fluidos se utiliza el método de elementos finitos de partículas (Particle Finite Element Method (PFEM)) mientras que los sólidos se solucionan con el método de elementos finitos (Finite Element Method (FEM)). Por lo tanto, se ramalla sólo las partes del dominio ocupadas por el fluido. Los campos de velocidad y presión se interpolan con funciones de forma lineales. Para poder analizar materiales incompresibles, la formulación ha sido estabilizada con una nueva versión del método Finite Calculus (FIC). La técnica de estabilización ha sido derivada para fluidos Newtonianos casi-incompresibles. En este trabajo, la estabilización con FIC se usa también para el análisis de sólidos hipoelásticos casi-incompresibles. En la tesis se dedica particular atención al estudio de flujo con superficie libre. En particular, se analiza en profundidad el tema de las pérdidas de masa y se muestra con varios ejemplos numéricos la capacidad del método de garantizar la conservación de masa en problemas de flujos en superficie libre. Además se estudia con detalle el condicionamiento del esquema numérico analizando particularmente el efecto del módulo de compresibilidad. Se presenta también una estrategia basada en el uso de un pseudo módulo de compresibilidad para mejorar el condicionamiento del problema. La formulación unificada ha sido validada comparando sus resultados numéricos con pruebas de laboratorio y resultados numéricos de otras formulaciones. En la mayoría de los ejemplos también se ha estudiado la convergencia del método. En la tesis también se describe una estrategia segregada para el acoplamiento de la formulación unificada con el problema de transmisión de calor. Además se presenta una simple estrategia para simular el cambio de fase. El esquema acoplado ha sido utilizado para resolver varios problemas de FSI donde se incluye la temperatura y su efecto. El esquema acoplado con el problema térmico ha sido utilizado con éxito para resolver un problema industrial. El objetivo del estudio era la simulación del daño y la fusión de la vasija de un reactor nuclear provocados por el contacto con un fluido altamente viscoso y a gran temperatura. En la tesis se describe con detalle el estudio numérico realizado para esta aplicación industrial.

keywords

1.3

pdftitle=Unified Lagrangian formulation for fluid and solid mechanics, fluid-structure interaction and coupled thermal problems using the PFEM pdfsubject= pdfauthor=Alessandro Franci pdfkeywords=#1

http://www.upc.edu Universitat Politècnica de Catalunya University Name

Doctoral Thesis

Unified Lagrangian formulation for fluid and solid mechanics, fluid-structure interaction and coupled thermal problems using the PFEM

Author:

[ Alessandro Franci]

Supervisors:

[ Prof. Eugenio Oñate Ibáñez de Navarra

Dr. Josep Maria Carbonell i Puigbó]

A thesis submitted in fulfilment of the requirements

for the degree of Doctor of Philosophy

in Structural Analysis

https://www.rmee.upc.edu Departament de Resistència de Materials i Estructures a l'Enginyeria

Barcelona, JanuaryFebruaryMarchAprilMayJune JulyAugustSeptemberOctoberNovemberDecember

dummy

1.3

–––––––––––––––––––––––––––––––––––––––Acknowledgements––––––––––––––––––––––––––––––––––––––– dummy toc I would like to thank all the people that in some way have helped me to complete this doctoral thesis. I would like to express my gratitude to Prof. Eugenio Oñate for many reasons. Thanks to him, four years ago I could come to Barcelona and start my PhD at CIMNE. He gave me the balanced doses of freedom and pressure and he encouraged me to test innovative technologies. I would like thank Prof. Oñate also for giving me the opportunity and the material support to spend a three-month period at the College of Engineering of Swansea. Special thanks go to Dr Josep Maria Carbonell. Without his support I could not realize this work. He devoted me innumerable hours of his time and he gave me technical and not technical suggestions that helped me to make the right choices during my PhD. Moltes gràcies Josep Maria! I would like to thank Prof. Javier Bonet for giving me the possibility to work with his group during my stay at Swansea. Special thanks go to Prof. Antonio Gil y Dr. Aurelio Arranz for their very kind treatment, for improving my knowledge of the Immersed Boundary Potential Method and for the useful discussions about FSI strategies. I really hope that this is just the first step for a future fruitful collaboration. I would like to thank also Mr. Kazuto Yamamura and the Nippon Steel & Sumimoto Metal Corporation for giving us the possibility to publish the PFEM simulation of the damages of a nuclear power plant pressure vessel. I also acknowledge the Agencia de Gestió d'Ajuts Universitaris i de Recerca (AGAUR) for the financial support. I would like to thank also Riccardo, Antonia, Jordi, Pablo and Lorenzo for their contributions for this thesis. This thesis is dedicated to my parents that still represent to me a fundamental and unique support. Grazie per tutti gli sforzi che avete fatto in questi anni e per le numerose visite a Barcellona. Questa tesi è per voi. Finally, I would like to thank Lucía simply for staying by my side during these years. Her smile has been the gasoline I used for dealing with all this work.

Contents –––––––––––––––––––––––––––––––––––––––Table of Contents––––––––––––––––––––––––––––––––––––––– dummy

List of Figures –––––––––––––––––––––––––––––––––––––––List of Figures––––––––––––––––––––––––––––––––––––––– dummy

List of Tables –––––––––––––––––––––––––––––––––––––––List of Tables––––––––––––––––––––––––––––––––––––––– dummy

1.5

Abbreviations –––––––––––––––––––––––––––––––––––––––Abbreviations––––––––––––––––––––––––––––––––––––––– dummy

Abbreviations

AL Augmented Lagrangian

ASGS Algebraic SubGrid Scale

FCM Fuel Containing Material

FEM Finite Element Method

FIC Finite Increment Calculus

FSI Fluid-Structure Interaction

GLS Galerkin Least-Squares

IBM Immersed Boundary Method

IFEM Immersed Finite Element Method

ISPM Immersed Structural Potential Method

LBB Ladyzenskaya-Babuzka-Brezzi

LFCM Lava-like Fuel Containing Material

NPP Nuclear Power Plant

OSS Orthogonal SubScale

PFEM Particle Finite Element Method

SPH Smooth Particle Hydrodynamics

SUPG Streamline Upwind Petrov Galerkin

TL Total Lagrangian

UL Updated Lagrangian

V Velocity

VMS Variational Multi-Scale

VP Velocity-Pressure

VPS Velocity-Pressure-Stabilized

1.3

–––––––––––––––––––––––––––––––––––––––Dedicatory–––––––––––––––––––––––––––––––––––––––

toc

``Dai diamanti non nasce niente,

dal letame nascono i fior."


``Nothing grows out of precious diamonds,

Out of dung the flowers do grow."

Fabrizio De Andrè

1 Introduction

Chapter 1. Introduction

1.1 Objectives

The objective of this work is to develop a unified formulation for the solution of the fluid and solid mechanics, fluid-structure interaction (FSI) and thermal coupled problems. The aim is to analyze the continuum in a unified manner trying to reduce at minimum the differences between the analysis of fluids and solids. For this purpose the numerical model has been planned in order to meet the specific requirements of solid and fluid mechanics and their approximation with the Finite Element Method (FEM) IEEEexample:Zinc1, but without limiting excessively the capability of the model. In fact the computational method should be capable to deal with critical problems such as those involving elasto-plastic solids, quasi-incompressible materials, free surface fluids and phase change due to heat transfer solutions.

According to these considerations, the computational model has been designed for using a stabilized Velocity-Pressure formulation. The formulation has been particularized for hypoelasto-plastic, compressible and quasi-incompressible solids and quasi-incompressible Newtonian fluids. The algorithm for the FSI problems has been inspired by the analogous unified strategy presented in IEEEexample:Idelshon. For the fluid phase, the Particle Finite Element Method (PFEM) IEEEexample:PFEM2006 has been used, while for the solid the classical FEM is adopted. The scheme has been stabilized with an updated version of the Finite Calculus (FIC) technique IEEEexample:nuestroPaper. For dealing with thermal coupled problems, this unified formulation has been extended to treat the heat transfer problem via a staggered scheme. For the modeling of the phase change, a simplified procedure has been developed. The whole formulation has been implemented in a C++ code.

1.2 State of the art

In this section, an overview of the numerical methods used for simulating a free surface fluid flow interacting with deformable solids is given. For the sake of clarity, the section is divided in three parts representing the main fields involved by this work. First the Eulerian and Lagrangian approaches for free surface fluid dynamics problems are presented. Then an overview of the stabilization for incompressible, or quasi-incompressible, material is given. Finally, the principal algorithms for solving FSI problems are described.

1.2.1 Eulerian and Lagrangian approaches for free surface flow analysis

Consider the description of the motion of a general continuum represented in Figure 1. The domain represents the body at the initial state at time while the domain represents the same body at time after deformation. The domain is called initial configuration, whereas is the current, or deformed, configuration. In order to describe the kinematics and the deformation of the body, the reference configuration has to be defined because the motion is defined with respect to this configuration.
Error creating thumbnail: File missing
Figure 1: Description of the motion of a general continuum body
In solid mechanics, the stresses generally depend on the history of deformation and the undeformed configuration must be specified in order to define the strains. Due to the history dependence, Lagrangian descriptions are prevalent in solid mechanics. In fluids however, the stresses do not depend on the history and it is often unnecessary to describe the motion with respect to a reference configuration. For this reason an Eulerian description represents the most reasonable choice. Furthermore, in problems where the fluid contours are fixed, Eulerian meshes are generally preferred to the Lagrangian ones. This is because Eulerian grids are fixed and they do not deform according to the fluid motion, as shown in Figure 2.
Error creating thumbnail: File missing
Figure 2: Motion description using an Eulerian mesh
Conversely, in the Lagrangian description, the mesh nodes coincide with the fluid particles and the discretization moves and deforms as the fluid flow (Figure 3).
Error creating thumbnail: File missing
Figure 3: Motion description using a Lagrangian mesh

Consequently, on the one hand the non-linear convective term disappears from the problem and on the other hand the mesh undergoes large distortions and it requires to be regenerated.

In the analysis of free surface flows, the detection of the free surface contours represents a crucial task. Its position is unknown and it has to be determined at each time increment in order to solve properly the boundary value problem. For these problems, the Lagrangian description may be preferred to the Eulerian one. In fact, with a Lagrangian approach the free surface is detected automatically by the position of the mesh nodes, while an Eulerian approach requires the implementation of a specific technology for this task.

Several strategies have been developed in the literature for tracking the free surface in an Eulerian framework. One of the earliest contributions was given by the so called marker and cell method IEEEexample:MarkerParticle. In this approach a set of marker particles that move according the flow are used to detect which regions are occupied by the fluid and which not. An evolution of this technique is the volume of fluid method IEEEexample:VOF. In this case the free surface boundary is detected using a scalar function that assumes the unit value in the fluid cells and the value zero in those ones with no fluid. The cells with an intermedium value are the ones that contain the free surface. Another possibility is the level set method IEEEexample:LevelSet. This technique is used in various fields, not only in continuum mechanics, and it allows for detecting shapes or surfaces on a fixed grid without making any parametrization of them. For this reason, this procedure has been also used for matching the free surface contour on an Eulerian mesh IEEEexample:antoniaLevelSet. A similar idea was used in IEEEexample:pabloPFEM2 where the position of the free surface is detected using a cloud of Lagrangian particles moving over an Eulerian mesh.

The free surface flows can be solved also using an hybrid Eulerian-Lagrangian technique. This is the so termed Arbitrary Lagrangian-Eulerian (ALE) approach IEEEexample:Zinc3. The aim of this method is to exploit the best features of the Eulerian and the Lagrangian procedures and to combine them. The mesh nodes can arbitrary be fixed or can move with the fluid IEEEexample:Donea. Generally, far from the moving boundaries a fixed Eulerian grid is used, while near to the interface the mesh moves according to the motion of the boundary IEEEexample:nargesALE. However if the boundary motion is large or unpredictable, also in the ALE methods the grid may suffer large distortions and may require a proper remeshing procedure IEEEexample:ALE.

In purely Lagrangian approaches, the mesh needs to be regenerated whenever a threshold limit for the distortion is exceeded. This is the basis of a particular class of Lagrangian finite element formulation called the Particle Finite Element Method (PFEM). The method was initially developed by the group of professors Idelsohn and Oñate IEEEexample:PFEM2004, IEEEexample:CIMNEPFEM. The PFEM treats the mesh nodes of the domain as particles which can freely move and even separate from the rest of the fluid domain representing, for instance, the effect of water drops. A mesh connects the nodes discretizing the domain where the governing equations are solved using a classical finite element method. These features make the PFEM the ideal numerical procedure to model and simulate free surface flows. In the last years, many scientific publications have shown the efficacy of the PFEM for solving free surface flow problems, see among others IEEEexample:PFEM2008,IEEEexample:cremonesiPFEM2,IEEEexample:PFEMtang. The PFEM has also been tested successfully in other kind of problems, such as fluid mechanics including thermal convection-diffusion IEEEexample:PFEMthermal,IEEEexample:PFEMmelting,IEEEexample:PFEMbedErosion2, multi-fluids IEEEexample:pabloMultiF, IEEEexample:PFEMmeschke, granular materials IEEEexample:PFEMzhangGranular, bed erosion IEEEexample:PFEMbederosion, FSI IEEEexample:PFEMfsi, IEEEexample:PFEMzhuFSI and excavation IEEEexample:PFEMexcavation.

Meshfree methods are other class of Lagrangian techniques. In this strategy the remeshing is not required because the governing equations are solved over a set of nodes without referring to a mesh. One of the first meshfree techniques is the Smooth Particle Hydrodynamics (SPH) method. This method was introduced independently by Gingold and Monaghan IEEEexample:SPH and Lucy IEEEexample:Lucy for the simulation of astrophysical problems such as fission of stars. SPH is a particle-based Lagrangian technique where discrete smoothed particles are used to compute approximate values of needed physical quantities and their spatial derivatives. The particles have assigned a characteristic distance, called ‘smoothing length’, over which their properties are ”smoothed” by a kernel function. A typical drawback of the SPH method is that it is hard to reproduce accurately the incompressibility of the materials. The SPH technique has been used successfully for solving fluid-structure interaction problems IEEEexample:SPHfsi.

1.2.2 Stabilization techniques

A FEM-based procedure may require to be stabilized when incompressible, or quasi- incompressible, materials are analyzed. For example from the finite element solution of the Navier Stokes equations numerical instabilities may arise from two sources. The first one is due to the presence of the convective term in the linear momentum equations. This term introduces a non-linearity in the equations and it needs a proper stabilization for solving high Reynolds number flows with the FEM IEEEexample:Donea. Furthermore, the orders of interpolation of pressure and velocity fields cannot be chosen freely but they have to satisfy the so called (LBB), or , condition IEEEexample:Brezzi. If the orders of interpolation of the unknown fields do not satisfy this restriction, a stabilization technique is required in order to avoid numerical instabilities, as the spurious oscillations of the pressure field.

It is well known that the weak form generated by Galerkin approximation leads to a less diffusive solution than the strong problem. So the main idea of many stabilization techniques consists of adding an artificial diffusion to the problem. The first attempt was made by Von Neumann and Richtmyer IEEEexample:VonNeumann. Their solution adds an artificial diffusion to the strong form of the problem. This technique introduces an excessive dissipation to the problem because the diffusivity is added in every direction. An important evolution of this approach was the Streamline Upwind Petrov Galerkin (SUPG) IEEEexample:SUPG. In this approach, the artificial diffusion is added by means of the test functions and only on the direction of the streamlines. Furthermore, this is performed in a consistent way: the stabilization terms vanish when the solution is reached. An extension of the SUPG method is the Galerking Least-Squares (GLS) method IEEEexample:GLS. The main difference between the two is that in the GLS method the stabilization terms are applied not only to the convective term but to all the terms of the equation. The Variational Multi-Scale (VMS) methods IEEEexample:VMM split the problem variables in a large-scale and a subscale terms. The large-scale terms represent the part of the solution that can be captured by the finite element mesh, while the subscale part consists of an approximate solution that has to be added to the large-scale term in order to obtain the correct solution. This idea represents the basis of other stabilization methods. Among these, the most largely used are the Algebraic Subgrid Scale Formulation (ASGS) and the Orthogonal Subscales (OSS) methods, respectively introduced in IEEEexample:ASGS and IEEEexample:OSS. Another efficient stabilization technique is the Finite Calculus (also termed Finite Increment Calculus) (FIC) approach (see among others IEEEexample:FIC1,IEEEexample:FICadt,IEEEexample:FICfsi2,IEEEexample:FICale ). This method has some analogies with the SUPG technique (for a comparison between these methods see IEEEexample:FICsupg). The FIC approach is based on expressing the equations of balance of mass and momentum in a space-time domain of finite size and retaining higher order terms in the Taylor series expansion typically used for expressing the change in the transported variables within the balance domain. In addition to the standard terms of infinitesimal theory, the FIC form of the momentum and mass balance equations contains derivatives of the classical differential equations in mechanics multiplied by characteristic distances in space and time. In this work an updated version of the FIC method has been derived and tested.

In solid mechanics a stabilization procedure may be required for the solution of problems involving incompressible, or quasi-incompressible solids. Situations of this type are common in forming processes or in the analysis of the rubber-type materials. Many of the stabilization procedures for fluid dynamics have been also used also for solid dynamics. For example, the VMS method has been applied in quasi-incompressible solid mechanics in IEEEexample:Chiumenti1, IEEEexample:ChiumentiInc, the OSS in IEEEexample:Chiumenti. In IEEEexample:GilPGsolid, IEEEexample:GilPGsolid2 a stabilized multi-field Petrov-Galerkin procedure is used. Finally, the application of the FIC in solid mechanics is reported IEEEexample:rojekFIC, IEEEexample:Onate.

1.2.3 Algorithms for FSI problems

Many approaches have been developed for solving FSI problems. Typically the computational techniques for FSI are distinguished in monolithic and staggered, or partitioned, approaches. In a monolithic approach the fluid and the solid domains are solved in a single system of equations (see among others IEEEexample:FSImono1, or IEEEexample:FSImono2). In this technique the flow of information between the solid and the fluid parts is implicitly performed by the procedure. On the contrary, in staggered schemes the fluid and the solid dynamics are solved separately and boundary conditions are transferred from a domain to the other at the interface. From the algebraic point of view, the solution is achieved by solving two different linear systems which are coupled by means of the boundary conditions defined along the interface. Within partitioned schemes, depending on the level of coupling between the fluid and the solid dynamics, a further classification can be done. In a weakly coupled segregated schemes when the transfer of information through the interface is performed only once for each time step, (for example see IEEEexample:FSIseg1), whereas in strongly coupled schemes this operation is performed within a convergence loop (as a reference see IEEEexample:FSIseg3). Clearly, a weakly coupled scheme has a lower computational cost but it can be used only when the interaction between the solid and the fluid domains is not strong or complex. Otherwise the algorithm may not find the correct numerical solution of the problem.

Partitioned methods allow the reutilization of existing solvers. Furthermore, the solid and the fluid solvers can be updated independently. Staggered schemes lead to smaller and better conditioned linear systems than the ones obtained with monolithic approaches. On the other hand, monolithic strategies are generally more stable and fast than staggered schemes, and they lead to a more accurate solution of FSI problems, because a stronger coupling is ensured.

Another general classification of the FSI algorithms is based upon the treatment of meshes. In conforming mesh methods, in order to allow the transfer of information, the fluid and the solid meshes must have in common the nodes along the interface. Consequently, if the position of the interface nodes changes in a domain, also the other domain must modify its grid in order to guarantee the conformity of the two meshes along the interface. On the contrary, in non-conforming mesh methods the interface and the related conditions are treated as constraints imposed on the model equations so that the fluid and solid equations can be solved independently from each other with their respective grids IEEEexample:FSInonC. This represents an important advantage because, typically, the mesh used for the fluid has an average size lower than the one used for the solid and so it is not necessary to refine the solid finite element grid near to the interface. However, non-conforming mesh algorithms are more complex to implement and it is not trivial to guarantee their robustness.

In the so-called (), the fluid is solved using an Eulerian grid and the solids are immersed on top of this mesh IEEEexample:PeskinIBM, IEEEexample:PeskinIBMuns. The interaction is ensured by penalizing the Navier-Stokes equations with the momentum forcing sources of the immersed structures. An evolution of the IBM is the () where the structure is modeled as a potential energy functional solved over a cloud of integration points that move within the fixed fluid mesh IEEEexample:GilISPM, IEEEexample:PeskinIBMuns. Also in the () IEEEexample:zhangIFEM, IEEEexample:GilIFEM the structure acts as a momentum forcing source for the fluid governing equations, but in this case a Lagrangian mesh is employed for the solid domains.

1.3 Numerical model

The aim of this work is to derive a finite element formulation capable to solve, through a unique set of equations and unknown variables, the mechanics of a general continuum. The term 'general continuum' refers to a domain that may include compressible and quasi-incompressible solids, fluids or both interacting together. For this reason, the formulation is termed . The Unified formulation is based on a mixed Velocity-Pressure Stabilized procedure and it has been implemented in a sequential code.

1.3.1 Reasons

There are many reasons for undertaking the above objective.

The first advantage of the Unified formulation is that it allows for solving fluid and solid dynamics problems by implementing and using a single calculation code.

Furthermore, if solids and fluids are solved via the same scheme, it is simpler to implement the solver for FSI problems because it is not required neither changing the variables, neither implementing the transfer of transmission conditions through the interface. With this formulation solids and fluids represent regions of the same continuum and they differ only by the specific values of the material parameters. As a consequence, the FSI solver requires a small computational effort and it can be implemented by introducing just a few specific functions. This will be explained in detail in the section dedicated to FSI problems.

Additionally, with the Unified formulation the most natural approach for solving FSI problem is the monolithic one. This brings in the further advantage that the coupling is ensured strongly and an iteration loop is not required, as for staggered procedures.

Finally, using the same set of unknowns for the fluid and the solid domains reduces the ill-conditioning of the FSI solver, because the solution system does not include variables of different units of measure.

1.3.2 Essential features

The Unified formulation is a compromise between a fluid and a solid formulation and it should be capable to satisfy the requirements of each problem mechanics and finite element approximations. In other words, in order to solve adequately fluids and solids with the same formulation, the solution procedure must take into account the constraints that both models impose.

Concerning the choice of the unknown variables, these are generally selected depending on the constitutive law of the materials. For Newtonian fluids, the Cauchy stress tensor is related directly to the deformation rate tensor. This implies that the velocities are the most useful unknowns in fluid mechanics. Conversely, in solid mechanics the displacements are generally preferred as variables, as the stresses are related to the deformations. However, for solids there also exist constitutive laws expressed in rate form, as for hypoelastic models.

For the analysis of incompressible, or quasi-incompressible, materials a mixed formulation is required in order to overcome the associated numerical instabilities, both in solids and fluids. In fluids, the most popular combination of unknown variables is the velocity-pressure scheme. In solid mechanics, there are several mixed approaches, as the displacement-stress IEEEexample:BrezziFortin, the displacement-strain IEEEexample:Lorenzo, the displacement-deformation gradient IEEEexample:GilPGsolid, the displacement-deformation gradient-jacobian IEEEexample:GilPGsolid2, the displacement-pressure IEEEexample:Chiumenti and the velocity-pressure formulation IEEEexample:PavelPaperA. The combination of unknowns is selected depending on the desired accuracy for the stress and strain fields and the associated computational cost. Among the mentioned mixed approaches, the displacement-pressure and the velocity-pressure formulations lead to the lowest computational cost, because the number of unknowns is smaller. However the velocity-pressure formulation has the further advantage that is the canonical scheme for fluid mechanics. Hence, according to these considerations, the Unified formulation has been based on a mixed velocity-pressure scheme.

Another key decision in the design of the Unified formulation concerns the choice of the description framework. Solids are typically solved using a (total or updated) Lagrangian description, while for fluids, due to the large deformations they undertake, the Eulerian framework is generally preferred. However, for free surface flows, as the ones treated in this work, a Lagrangian description has the important advantage that the free surface boundaries are detected automatically, because the particles of the fluid coincide with the nodes of the mesh. The price for this is that a remeshing procedure is required in order to avoid the excessive distortion of the mesh. According to these considerations, the Unified formulation has been implemented using an updated Lagrangian description.

In this work, the Particle Finite Element Method (PFEM) IEEEexample:PFEM2006 has been used for solving the fluid domain only. The PFEM is a Langrangian numerical technique that treats the mesh nodes as moving points which can freely move and even separate from the domain to which they are initially attached representing, for instance, the effect of water drops. The PFEM is based on a remeshing procedure that efficiently combines the Delaunay tessellation and the Alpha Shape method IEEEexample:Edelsbrunner. In order to reduce the computational cost associated to the remeshing step, all the unknowns are stored in the nodes and linear shape functions are used for the finite element interpolation. Conversely, the solid domain keeps a fixed mesh during the whole analysis. The reason for this is that in non-linear solid mechanics it is required to preserve not only the nodal unknowns but also elemental information as the stress state of the previous time step. A remeshing implies the elimination of the previous discretization and the creation of new elements. In order to keep the elemental information in the remeshing step, a projection procedure from the elements of the previous mesh to its nodes, and from these to the elements of the new mesh is required. These operations may introduce an interpolation error in the scheme that affects the solution of solid mechanics problems. In the analysis of Newtonian fluids all the information can be stored in the nodes, so the remeshing does not affect the accuracy of the scheme. For these reasons, in this work the PFEM is used only for the fluid, while solids are solved via the classical FEM.

In order to improve the conditioning of the solution system, a Gauss-Seidel segregated scheme is used. This means that the problem is solved iteratively and separating the unkwown variables, the velocities and the pressure. The solution method consists of two steps. In the first one the linear momentum equations are solved for the velocity increments and considering the pressure of the previous non-linear iteration in the residual term. In the second step, the continuity equation is solved for the pressure in the updated configuration using the velocities computed at the first step. This two-step algorithm leads to a smaller and better conditioned linear system than using a monolithic approach. A crucial point of this scheme is the derivation of the tangent matrix of the linear momentum equations. The velocity increments are solved by condensing the pressure. According to this procedure, the tangent matrix for the linear momentum equations of the mixed velocity-pressure formulation does not differ from the one of a velocity formulation. Furthermore, the derivation of the tangent matrix for the velocity formulation is easier, as it does not involve the pressure. For these reasons, the first step towards the unified scheme is the derivation of the Velocity formulation. The mixed Velocity-Pressure formulation will be derived by exploiting the linearization performed for the Velocity formulation. In IEEEexample:LagrangianTangent another procedure for the derivation of the exact tangent matrix for an Updated Lagrangian scheme is described.

The condensation of the pressure in the tangent matrix of the linear momentum equations may induce the ill-conditioning of the linear system for the analysis of quasi-incompressible materials. This ill-conditioning emanates from the volumetric counterpart of the tangent matrix that is governed by the bulk modulus, that typically has high values that may compromise the conditioning of the linear system. In this work, a thorough study of the conditioning drawbacks associated to this scheme has been performed and a useful and easy to implement technique to overcome these inconveniences has been tested and validated.

The numerical method proposed in this work is an improvement versus other numerical schemes for quasi-incompressible materials where an arbitrarily defined pseudo-bulk modulus was used in both the linear momentum and the mass conservation equations IEEEexample:PavelPaperA, IEEEexample:PavelPaper, IEEEexample:PavelPaperB.

The method can be also related to the so-called Augmented Lagrangian (AL) procedures for solving the Navier-Stokes equations for incompressible IEEEexample:Benzi1,IEEEexample:Benzi2,IEEEexample:Elman,IEEEexample:Fortin,IEEEexample:Vincent and weakly compressible IEEEexample:Bresch,IEEEexample:Vinay flows.

In order to deal with incompressible (or quasi-incompressible) materials a mixed formulation is required. In this work both the velocities and the pressure are interpolated using linear shape functions. This combination does not fulfil the condition IEEEexample:Brezzi and the problem needs to be stabilized. The required stabilization is ensured using an updated version of the FIC technique IEEEexample:nuestroPaper applied to the mixed Velocity-Pressure formulation. The FIC stabilization method has a small intrusivity because its terms only affect the continuity equation. In fact for solving the linear momentum equations the same scheme as for the (not stabilized) mixed Velocity-Pressure formulation is used. It will be shown that the PFEM-FIC stabilized procedure guarantees a good accuracy for the mass conservation of the free surface flows. The stabilization procedure is derived for quasi-incompressible fluids, but it can be easily extended to quasi-incompressible hypoelastic solids.

The derivation of the Unified stabilized formulation is carried on trying to maintain as much as possible the generality of the scheme. For this reason the constitutive relations are introduced only as the last step. The Velocity (V) and the mixed Velocity-Pressure (VP) formulations are derived first for a general material and then particularized for specific constitutive laws. For quasi-incompressible materials only the mixed Velocity-Pressure Stabilized (VPS) formulation can be used while for compressible materials the Velocity and the standard (not stabilized) mixed Velocity-Pressure formulations are both suitable. In particular, for compressible solids, the hypoelastic model has been chosen.

The FSI problem is solved in a monolithic way, hence solids and fluids are solved by the same linear system. For the solid domain it is possible to choose which formulation to use. Depending on the problem, one may chose the V, the VP or the VPS element.

The Unified formulation for FSI problems can be easily coupled with the heat transfer problem in order to solve coupled thermal mechanical problem. The coupling is performed via a staggered scheme. The effect of the heat is taken into account by considering the material properties depending on the temperature and including in the strain tensor the deformations induced by the temperature. With the PFEM Unified formulation also phase change problems can be modeled.

The Unified stabilized VP formulation with thermal coupling has been used for the analysis of an industrial application. The study concerned the analysis of the damages caused by the dropping of a volume of corium at high temperature on the structure of the pressure vessel in a Nuclear Power Plant (NPP).

1.3.3 Outline

This text is split in the following chapters.

In the next chapter two velocity-based FEMs are derived for a general compressible material. In the first section the Velocity formulation is derived and the incremental solution scheme is given. Then, the standard mixed Velocity-Pressure formulation is obtained as an extension of the Velocity formulation. The constitutive laws are not introduced for any scheme. In the third section the formulations are adapted for hypoelastic-plastic compressible solids. The chapter ends with several validation examples for non-linear solid mechanics.

In Chapter 3, the FIC stabilization strategy is introduced in the mixed stabilized VP scheme. This scheme is adapted for quasi-incompressible Newtonian fluids and hypoelastic solids, in this order. Then the free surface problem is studied in detail. First the PFEM is described highlighting the advantages and the disadvantages of the method. A useful technique for modeling the slip conditions in Lagrangian flows is also explained. Next the mass conservation feature and the conditioning of the scheme are analyzed in detail. Concerning the former point, it will be shown that the PFEM-FIC formulation guarantees a good preservation of the mass in free surface flows. Regarding the latter point, a practical strategy for improving the conditioning of the linear system is explained and tested. At the end of the chapter several validation examples for quasi-incompressible Newtoninan fluids and hypoelastic solids are given.

The FSI solution strategy is described in Chapter 4. The monolithic scheme for coupling the mechanics of fluids and solids is explained in detail. Finally several validation examples of FSI problems are given.

Chapter 5 is dedicated to couple the Unified formulation for FSI with the heat problem. In the first section the heat problem is introduced and discretized using the FEM. Then the coupling strategy is explained and validated with numerical examples. Next, the procedure for modeling phase change problems is described and an explicative numerical example is presented.

Chapter 6 is fully devoted to the industrial application of the Unified stabilized and thermally coupled strategy. The damages of the pressure vessel structure caused by the fall of a volume corium at high temperature are studied using two simplified models where the solid structure melts due to the heat transfer from the corium.

In Chapter 7 the innovative contributions of the work are summarized and the lines of research opened by this thesis are presented.

1.4 Publications

The following publications have emanated from the work carried out in this doctoral thesis

2 Velocity-based formulations for compressible materials

Chapter 2. Velocity-based formulations for compressible materials

In this chapter two velocity-based finite element formulations for compressible materials are presented, namely the Velocity (V) and the mixed Velocity-Pressure (VP) formulations. For both schemes the linear momentum equations are solved iteratively for the velocity increments. The linearization of the governing equations is performed without specifying any constitutive law. The aim of this chapter is to maintain as much as possible the generality of the algorithms, leaving the formulations open to different material models. It will be shown that the only requirement demanded to the constitutive laws is that the rate of stress must be linearly related with the rate of deformation.

There are several reasons that justify the presentation of the Velocity formulation. First of all, the tangent matrix of the linear momentum equations for the Velocity formulation holds also for the mixed Velocity-Pressure formulation but the linearization procedure is easier because the pressure does not appear. Furthermore, the Velocity formulation is useful for making some interesting and didactic comparisons with the mixed formulations.

After deriving the solution scheme for the Velocity formulation, the mixed Velocity-Pressure method is presented. The governing equations are the linear momentum and the linear pressure-deformation rate equations. The latter is called continuity, or mass balance, equation for its similarity with the incompressibility constraint of the Navier-Stokes problem. As for the velocity scheme, in the mixed formulation the constitutive law is not specified. In Section 2.3 the hypoelastic-plastic model is presented and this constitutive law is inserted in both the Velocity and the mixed Velocity-Pressure schemes. The incremental solution scheme is explained in detail for both formulations. At the end of the chapter some validation examples for hypoelastic-plastic solids in statics as in dynamics are given.

2.1 Velocity formulation

In this section, the velocity formulation for solving transient problems for a general continuum is derived. The governing equations are the linear momentum equations and they are derived in the updated Lagrangian (UL) framework. This means that the governing equations are integrated over the unknown configuration (the so-called updated configuration). As a consequence, the space derivatives for the UL description are computed with respect to the spatial coordinates.

2.1.1 From the local form to the spatial semi-discretization

In this section the spatial semi-discretization of the linear momentum equations is derived.

For a general continuum, the local form of the linear momentum equations using the UL description reads

(1)

where is the density of the material, are the velocity vector, is the Cauchy stress tensor and is the body force vector. The variables within the brackets are the independent variables: X are the Lagrangian or material coordinates vector, x the Eulerian or spatial coordinates vector and is the time. For simplicity, in what follows the independent variables are not specified. The spatial and material coordinates are related through the motion tensor as

(2)

The set of governing equations is completed by the following conditions at the Dirichlet () and Neumann () boundaries

(3)

(4)

where and are the prescribed velocities and the prescribed tractions, respectively, and is the unit normal vector.

In the following summation of terms for repeated indices is assumed, unless otherwise specified.

The spaces for the trial and test functions are defined, respectively, as

(5)

(6)

Multiplying Eqs.(1) by the test functions and integrating over the updated configuration domain, the following global integral form is obtained

(7)

where the symbol represents the material time derivative.

Integrating by parts the term involving in Eq.(7) and using the Neumann boundary conditions (4) yields the weak variational form of the momentum equations as

(8)

Eq.(8) is the standard form of the Principle of Virtual Power [IEEEexample:BelyBook].

The spatial discretization is introduced using the classical FEM-Galerkin procedure [IEEEexample:Zinc1]. Hence both the trial and the test functions are interpolated in space in terms of their nodal values by means of the same shape functions

(9)

where, assuming the use of simplicial elements, for 2D/3D problems is the number of the nodes of the element, denotes a nodal value, the capital subscript specifies the node and the lower case subscript represents the cartesian direction. In this work, linear shape functions have been used for .

Since Eq.(9) must hold for all the test functions in the interpolation space, introducing the spatial discretization (9) into Eq.(8), the spatial semi-discretized form of the momentum equations in the UL framework for the node reads

(10)

where , and are the dynamic, internal and external forces, respectively, expressed in the UL framework.

For convenience, the semi-discretized form of the momentum equations in the total Lagrangian (TL) framework is also presented here. This is written as [IEEEexample:BelyBook]

(11)

where is the first Piola-Kirchhoff stress tensor, or the nominal stress tensor, and , and are the dynamic, internal and external forces, respectively, expressed in the TL framework. All the variables with vectors subscript refer to the last known configuration. Note that Eq.(11) can be obtained from Eq.(10) by pull back transformations on all its terms [IEEEexample:BelyBook].

For the sake of clarity in the notation, the terms referred to the TL description are denoted with the left index . Unless otherwise specified, the variables belong to the UL description.

2.1.2 Time integration

In this work, the kinematic variables have been integrated in time using a second order scheme. In particular, the implicit Newmark's integration rule has been adopted. For the general case, this rule states that accelerations and displacements are computed as

(12)

Where and are the so-called Newmark's parameters [IEEEexample:BelyBook]. This time integration scheme is unconditionally stable if the following relation holds

(13)

In the present work, the Newmark's parameters chosen are and .

Replacing the numerical values of the constants in Eq.(12) yields

(14)

(15)

2.1.3 Linearization

Although the problem is set out in the UL framework, the linearization for the velocities increment of the momentum equations is performed first on the TL semi-discretized form (11). The UL linearized form is obtained by performing a push-forward transformation on the TL form. This is justified by the easier derivation of the tangent matrix in the TL framework. In fact, in Eq.(11) the only variable that depends on time is the nominal stress , while in the UL form (10) the time-dependent variables are the updated domain , the Cauchy stress tensor and the spatial derivatives . For the sake of clarity, the linearization of the internal and dynamic forces will be performed separately.

From Eq.(11) the internal forces in the TL description are defined as

(16)

The constitutive relation is expressed in rate form. Hence it is more convenient to perform the linearization of the material derivative of the internal forces and then integrate for the time step increment . The material time derivative of (16) is

(17)

The infinitesimal increment of Eq.(17) is

(18)

The first Piola-Kirchhoff stress tensor is not typically used because it is not symmetric and its rate is a non-objective measure. For these reasons, in the TL framework it is more convenient to work with the second Piola-Kirchhoff stress tensor and its rate. These stress rate measures are related each other through the following relation

(19)

where is the deformation gradient tensor defined as

(20)

Substituting Eq.(19) into (18) yields

(21)

Eq.(21) shows that the increment of the material time derivative of the internal forces can be split into material and geometric parts, and , respectively. The former accounts for the material response through the rate of the second Piola-Kirchhoff stress tensor and the second term is the initial stress term that contains the information of the updated stress field. Note that, up to now, no constitutive relationships have been introduced and the above derivation holds for a general continuum.

From Eq.(21), the material part of the material time derivative of the internal forces reads

(22)

For the derivation of the material tangent matrix, the constitutive relation between the stress and the strain measures is required. In order to maintain the formulation as general as possible, the stress rate measure is related to the deformation rate through a tangent constitutive tensor, as

(23)

where is a fourth-order tensor and is the Green-Lagrange strain tensor. Eq.(23) can also be expressed in Voigt notation as

(24)

where denotes a vector with components . As it will be shown in the following sections, Eq.(23) can represent both a Kirchhoff solid material and a Newtonian fluid. If different constitutive laws are used, Eq.(23) should be modified accordingly in order to derive the material part of the tangent matrix.

The Green-Lagrange strain tensor can be expressed in terms of the nodal velocities as

(25)

In Voigt notation

(26)

where for a plane strain problem

(30)

Substituting Eq.(25) in Eq.(23), yields

(31)

Similarly, substituting Eq.(26) in Eq.(24)

(32)

Substituting Eq.(31) into Eq.(22) yields

(33)

where is the -component of the velocity of node . In Voigt notation Eq.(33) reads

(34)

In order to obtain the increment of the internal forces, the material time derivative of the internal forces increment is integrated over a time step increment as

(35)

From Eq.(35) and Eq.(33), yields

(36)

and, similarly, from Eq.(35) and Eq.(34)

(37)

From Eq.(36) and Eq.(37), the material tangent matrix in the TL description can be computed as

(38)

(39)

The material tangent matrix for the UL framework is obtained by applying a push-forward transformation on each term of Eq.(38) and integrating over the updated domain . The following relations hold

(40)

(41)

(42)

where is the tangent moduli corresponding to the material time derivative of the Kirchhoff stress tensor and is the tangent moduli for the rate of the Cauchy stress . The rate of the Cauchy stress tensor is related to the rate of deformation through the fourth-order tensor by the following expression

(43)

Substituting Eqs.(40-42) into (38) and using the minor symmetries, the material tangent matrix for the UL reads

(44)

Using Voigt notation, the same matrix reads

(45)

For the node of a 3D element, matrix is

(51)

The geometric tangent matrix for the UL framework is derived next using the same procedure adopted for the material components. Hence, first the linearization is performed using the TL form and then the UL tangent matrix is obtained by performing the required transformation over the TL terms.

From Eq.(21)

(52)

where the rate of the deformation gradient is defined as

(53)

Substituting Eq.(53) into Eq.(52), the geometric components of the internal power in the TL description can be written as

(54)

Integrating Eq.(54) on time for a time step increment yields

(55)

From Eq.(55) the geometric tangent matrix is obtained as

(56)

In order to recover the UL form, the Piola identity has to be recalled,

(57)

The geometric tangent matrix in the UL framework is obtained by substituting Eqs.(40), (41) and (57) into (56) and using the symmetries. This leads to

(58)

or also

(59)

where is the second order identity tensor and matrix for 3D problems is

(60)

The dynamic component of the tangent matrix in the UL description can be derived directly from the UL dynamic term of Eq.(10). This reads

(61)

Eq.(61) has to be discretized in time with the purpose of replacing the accelerations by the velocities using the time integration scheme described in Eq.(14).

Introducing Eq.(14) into (61) and differentiating for the increment of velocities, the dynamic components of the tangent matrix are obtained as

(62)

Or also

(63)

2.1.4 Incremental solution scheme

The problem is solved through an implicit iterative scheme. At each iteration the velocity increments are computed as

(64)

where is the tangent matrix computed as the sum of the internal, the geometric and the dynamic components, given by Eqs. (45), (59) and (63) respectively, as

(65)

is the residual and it is computed from Eq.(10) as

(66)

In Eq.(66) is the residual of the momentum equations referred to node and the cartesian direction . Note that the Cauchy stress tensor still appears in its 'continuum' form because up to now it has not been written as a function of the nodal unknowns. This is done for keeping the generality of the formulation. Only after the introduction of the constitutive laws, it will be possible to compute the Cauchy stress tensor as a function of the nodal unknowns.

In Box 1 the iterative solution incremental scheme of the velocity formulation for a generic time interval of duration is described.

2.2 Mixed velocity-pressure formulation

In this work, the mixed Velocity-Pressure formulation is derived as an extension of the Velocity formulation presented in the previous section. The governing equations are the linear momentum equations and the linear relation between the time variation of pressure and the volumetric strain rate. The latter represents a limit for the generality of this formulation because only materials which constitutive law satisfies this relation can be modeled through the mixed Velocity-Pressure formulation.

The problem is solved using a two-step Gauss-Seidel partitioned iterative scheme. First the momentum equations are solved in terms of velocity increments and including the (known) pressures at the previous iteration in the residual expression. Then the continuity equation is solved for the pressure using the updated velocities computed from the momentum equations. It will be shown that using this not intrusive scheme, it is possible to take advantage of most of the velocity formulation derived in the previous section. In particular, the incremental velocity scheme for the momentum equations (Box 1) and the structure of the tangent matrix (65) hold also for the mixed Velocity-Pressure formulation. The same linear interpolation has been used for the velocity and the pressure fields. It is well known that, for incompressible (or quasi-incompressible) problems, this combination does not fulfill the condition [IEEEexample:Brezzi] and a stabilization method is required. However, as mentioned in the section devoted to the velocity formulation, the aim of this part is to keep the formulation as general as possible without referring to a specific material. Hence only in the next chapter, when the mixed Velocity-Pressure formulation is used for solving quasi-incompressible problems, the required stabilization will be introduced in the scheme.

2.2.1 Quasi-incompressible form of the continuity equation

Mixed formulations are often used for dealing with incompressible materials. In these problems it is useful to write the stress and the strain measures as the sum of deviatoric and hydrostatic, or volumetric, parts. Hence the Cauchy stress tensor is decomposed as

(67)

with

(68)

where and are the deviatoric and the hydrostatic parts of the Cauchy stress tensor, respectively. The pressure is defined positive in the tensile state and equal to the hydrostatic parts of the Cauchy stress tensor . Hence

(69)

The Cauchy stress tensor can be computed as

(70)

The same decomposition is done for the spatial strain rate tensor . So

(71)

with

(72)

where and are the deviatoric and the hydrostatic parts of the strain rate tensor, respectively. The strain rate tensor is computed from the velocities as

(73)

The volumetric strain rate is defined from Eqs.(68) and (73) as

(74)

The closure equation for the mixed Velocity-Pressure formulation is the linear relation between the change in time of the pressure and the volumetric strain rate. This reads as

(75)

where is a parameter that depends on the constitutive equation. Typically is the bulk modulus of the material.

In conclusion, the local form of the whole problem for the mixed Velocity-Pressure formulation is formed by the linear momentum equations (Eq.(1)) and the pressure-strain rate relation given by Eq.(75). The linear momentum equations have been already discretized and linearized for the increments of velocities in the previous section devoted to the Velocity formulation. That form holds also for the mixed formulation. So, in the following, only the discretization of Eq.(75) is given. This equation is a restriction on the generality of the constitutive laws that can be analyzed with the mixed Velocity-Pressure formulation. It will be shown that the constitutive models for hypoelastic solids and quasi-incompressible Newtonian fluids fulfill this relation. In fluid dynamics, Eq.(75) represents the , or , equation for quasi-incompressible fluids. In fact, Eq.(75) with is the canonical form of the continuity equation of the Navier-Stokes problem. For this reason, from here on Eq.(75) will be called 'continuity equation'.

Multiplying Eq.(75) by arbitrary test functions (with dimensions of pressure), integrating over the analysis domain and bringing all the terms at the left hand side gives

(76)

Both the trial and the test functions for the pressure are interpolated in space using the same shape functions .

(77)

where for 2D/3D problems is the number of the nodes of the simplex. In this work, linear shape functions have been used for , as for the velocity.

Combining Eq.(77) with Eq.(76) and solving for all the admissible test functions q, yields

(78)

Regarding the time integration a first order scheme has been adopted for the pressure. Hence, for a time interval of duration the first and the second variations in time of the pressure are computed as

(79)

(80)

Introducing Eq.(79) in Eq.(78), the discretized form of the continuity equation is

(81)

where

(82)
(83)

(84)

where has been defined in Eq.(51) and .

2.2.2 Solution method

In the mixed Velocity-Pressure formulation the problem is solved through a partitioned iterative scheme. Specifically, the linear momentum equations are solved for the velocity increments as in the Velocity formulation. On the other hand, the continuity equation is solved for the pressure in the updated configuration using the velocity field computed with the linear momentum equations. In order to guarantee the coupling between the continuity equation and the linear momentum equations (or equally between the pressure and the velocities) the pressure must appear in the right hand side of the linear momentum equations. For this purpose the Cauchy stress tensor must be computed as the sum of its deviatoric part and the pressure, as Eq.(70). Otherwise, the Velocity-Pressure formulation would be uncoupled and totally equivalent to the Velocity formulation.

In conclusion, for a general time interval of duration the following linear systems are solved for each iteration

(85)

(86)

where is the same tangent matrix as for the Velocity formulation (Eq.(65)) and the residual is computed using the pressure of the previous iteration and the deviatoric part of the Cauchy stress as

(87)

In Box 2, the iterative incremental solution scheme for a generic continuum via the mixed velocity-pressure formulation is shown for a time interval .

2.3 Hypoelasticity

Using the definition of Truesdell [IEEEexample:Truesdell2], a hypoelastic body is a material which may soften or harden in strain but in general has neither preferred state nor preferred stress. The hypoelastic laws were created with the purpose of transferring the linear theory of elasticity from the small to the finite strains regime [IEEEexample:Truesdell2]. In [IEEEexample:Truesdell3] a deep dissertation about the differences between elasticity and hypoelasticity is given.

A hypoelastic body is defined by the constitutive equation [IEEEexample:Truesdell1]

(88)

In the rate theory it is crucial to guarantee the and the , or , of the scheme. A material is frame invariant if its properties do not depend on the change of observer. An objective constitutive equation is defined to be invariant for all changes of the observer [IEEEexample:Holzapfel]. For guaranteeing the frame indifference, the constitutive law has to be isotropic [IEEEexample:TruesdellLibro]. This represents a constraint for hypoelastic models. This limitation is even more severe if also plasticity is included in the model. In fact, for hypoelastic-plastic materials, also the yield condition is required to be isotropic [IEEEexample:SimoHughes].

The stress rate cannot be computed simply as a material derivative because it leads to a non-objective measure of stress [IEEEexample:BelyBook]. In particular, rigid rotations may originate a wrong state of stress if the stress rate is computed as the material time derivative of the Cauchy stress [IEEEexample:BelyBook]. However, many objective measures of rate of stress are available. The most common ones are the Truesdell's and Jaumann's Cauchy stress rate measures. From here on, an objective rate measure will be denoted by the upper inverse triangle index .

Most of hypoelastic laws relate linearly the stress rate to the rate of deformation. Hence, Eq.(88) is now rewritten in the following form

(89)

where is the Cauchy stress rate tensor, is the tangent moduli tensor and is the deformation rate tensor.

From Eq.(89) it can be deduced that this class of hypoelastic materials has a rate-independent and incrementally linear and reversible behavior. So, as for the elastic materials, in the small deformation regime, the strains and the stresses are totally recovered upon the unloading process. Nevertheless, for large deformations the hypoelastic laws do not guarantee that the work done in a closed deformation path is zero [IEEEexample:Prager]. However this error can be considered negligible if the elastic deformations are small with respect to the total deformations [IEEEexample:BelyBook]. For this reason the hypoelastic laws are often used for describing the elastic part of elastic-plastic materials: the plastic deformations in fact represent usually the largest part of the overall deformations.

The Jaumann measure for the rate of the Kirchhoff stress tensor is

(90)

where the tangent moduli fourth-order tensor for the Jaumann measure of the Kirchhoff stress rate is

(91)

where and are the Lamé constants and they are computed from the Young modulus and the Poisson ratio as

(92)

(93)

and is the fouth-order symmetric identity tensor defined as

(94)

Separating the volumetric from the deviatoric part, yields

(95)

where is the bulk modulus and it is computed from the Lamé parameters as

(96)

and is the fouth-order tensor computed as

(97)

The Jaumann stress rate measure is defined as

(98)

where is the Jaumann's tangent moduli tensor.

For a anisotropic material the tangent moduli for the Jaumann rate depends on the stress state and it is related to as follows

(99)

Instead, for isotropic materials, the Jaumann's tangent moduli tensors for the Cauchy stress rate and for the Kirchhoff stress rate are identical [IEEEexample:BelyBook]. So can be computed as

(100)

or equally

(101)

For a 3D problem tensor is

or, equally,

A material is said to be isotropic if its behavior is uniform in all directions, so it has no preferred orientations or directions. Many materials, such as metals and ceramics, can be modeled as isotropic for small strains [IEEEexample:BelyBook]. From the computational point of view, an isotropic constitutive law is much easier to manage than an anisotropic one and it has a lower computational cost. Isotropic laws are preferred, for example, for their symmetry properties. In fact the anisotropic tangent moduli (Eq.(99)) is not symmetric while, the isotropic one (Eq.(101)) has both minor and major symmetries, in fact

(102)

(103)

For all these reasons, in this work the isotropic law has been used for the hypoelastic model.

The tangent moduli will be introduced into the material part of the global tangent matrix Eq.(45) for the Velocity and the mixed Velocity-Pressure formulations indifferently.

As it has already pointed out, the Cauchy stress rate does not coincide with the material derivative of the Cauchy stress tensor. In fact, the following relation holds between both measures

(104)

where is a tensor that accounts for the rotations and it is defined as

(105)

where is the spin tensor defined as

(106)

In this work the tensor is computed at the end of each time step. Discretizing in time Eq.(104) for the time step interval and expanding the Cauchy stress rate, yields

(107)

can be viewed as a correction of the Cauchy stress tensor. So it can be related to the Cauchy stress tensor of the previous time step as follows

(108)

Replacing Eq.(108) in Eq.(107), yields

(109)

Substituting in Eq.(109) the relation for using Eq.(101), yields

(110)

Hence,

(111)

The first and the second terms of the right hand side of Eq.(111) represent the increment in the time step of the volumetric and deviatoric parts of the Cauchy stress tensor. From Eq.(111) it can be deduced that for isotropic hypoelastic solids described using the Jaumann measure, the following relation holds

(112)

Eq.(112) will be used as the closure equation of the mixed Velocity-Pressure formulation for hypoelastic solids. Note that Eq.(112) has the same structure as Eq.(75) analyzed in the previous section. From Eq.(112) using linear shape functions and integrating on time the pressure with a first order scheme, the same matrix form of Eq.(81) obtained for a general material is obtained.

In conclusion the updated stresses can be computed using the velocities only or both the pressure and the velocities, as follows

(113)

(114)

Eqs.(113) and (114) will be used for computing the Cauchy stress tensor in the Velocity and mixed Velocity-Pressure formulations, respectively.

2.3.1 Velocity formulation for hypoelastic solids

The solution scheme for hypoelastic solids is the one derived in Section 2.1.4 for a general continuum. The only modifications required are the definition of the tangent moduli in matrix (Eq.(45)) and the computation of the Cauchy stress tensor from the nodal velocities according to the hypoelastic model. This tensor appears in the geometric part of the tangent matrix (Eq.(59)) and into the residual (Eq.(66)). The tangent moduli tensor is taken from the Jaumann isotropic description and it is the tensor of Eq.(101). Concerning the Cauchy stresses, these are computed with Eq.(113). The finite element implemented with this hypoelastic velocity formulation is named V-element.

The iterative solution incremental solution scheme for hypoelastic solids using the velocity formulation for a generic time interval is given in Box 3.

All the matrices and vectors in Box 3 are collected in Box 4

2.3.2 Mixed Velocity-Pressure formulation for hypoelastic solids

The solution scheme is like the one presented in Section 2.2.2. As already explained, the tangent matrix of the mixed Velocity-Pressure formulation is the same as for the Velocity formulation. However the Cauchy stress tensor is computed from the nodal velocities and the nodal pressures using Eq.(114). The governing equations are Eqs.(85-86) particularized with the material parameters of a hypoelastic solid. The finite element implemented with this hypoelastic mixed velocity-pressure formulation is called VP-element.

In Box 5, the iterative solution incremental scheme for hypoelastic solids using the mixed velocity-pressure formulation is given for a generic time interval .

All the matrices and vectors that appear in Box 5 are collected in Box 6

Note that for the mixed velocity-pressure formulation the material part of the tangent matrix is defined by the same tangent moduli of the velocity scheme ().

In the mixed formulation, the momentum and the continuity equations can be easily decoupled. This is obtained by computing the Cauchy stress tensor using the velocities only (Eq.(113)) and not as the sum of its deviatoric part and the pressure (Eq.(114)). The uncoupled mixed velocity-pressure formulation is totally equivalent to the velocity formulation. In fact, although the pressures are still computed by solving the continuity equation, they do not appear in the momentum equations and, hence, they do not affect the solution for each the time step.

2.3.3 Theory of plasticity

The theory of plasticity is dedicated to those solids that, after being subjected to a loading process, may sustain permanent () deformations when completely unloaded [IEEEexample:Peric]. The plasticity is defined if the permanent deformations of the material do not depend on the rate of application of the loads. The materials whose behavior can be adequately described by this theory are called materials.

Elastic-plastic laws are characterized for being path-dependent and dissipative. The stresses cannot be computed as a single-valued function of the strains because they depend on the entire history of the deformation [IEEEexample:BelyBook].

The most important properties of the theory of plasticity can be summarized as follows:

  1. The increments of strain are assumed to be additively decomposed into an elastic (reversible) part and a plastic (irreversible) part , such that
    (115)
  2. There exists an where the behavior of the material is purely elastic and permanent deformations are not produced;
  3. The delimits the elastic domain. It governs the onset and the continuity of the plastic deformations and it is a functions of the state of stress and of the internal variables . So
    (116)
  4. The yield function cannot be positive: it is negative when the stress state is below the yield value and null otherwise (: );

  5. The plastic strain increments are governed by the so called ;
  6. is the and it is positive for a plastic loading and equal to zero for elastic loading or unloading;
  7. The loading-unloading process is described by the Khun-Tucker conditions
    (117)
  8. The third condition can also be expressed in the rate form through the so-called , . For plastic loading () the stress state lies on the yield surface (), instead for elastic loading or unloading the yield condition is not reached () and there is not plastic flow ().

2.3.3.1 Hypoelastic-plastic materials

Hypoelastic-plastic models are typically used when the elastic strains represent only a small part of the total strains. In other words, when the plastic strains are much larger than the elastic ones. This is because of the inaccuracy of the hypoelastic models in the large strain regime. However, if the elastic strains are small, the energy error introduced by the hypoelastic description of the elastic response is limited and can be considered adequate [IEEEexample:BelyBook].

Depending on the problem, the yield function can be based on a different constitutive model. For example, for soil plasticity the Drucker-Prager model is the most used, while for porous plastic solids the Gurson model is more adequate. In this work, the flow model is used. This model is particularly indicated for the metal plasticity [IEEEexample:BelyBook].

The hypoelastic-plastic model described in this section has been taken from [IEEEexample:BelyBook]. According to the von Mises criterion [IEEEexample:VonMises], plastic yielding begins when the stress deviator invariant reaches a critical value. The stress deviator invariant is defined as

(118)

The key assumption of the von Mises model is that the plastic flow is not affected by the pressure but only by the deviatoric stress. This hypothesis has been experimentally verified for metals [IEEEexample:Bridgman]. For this reason the von Mises model is called to be .

A yield function for the von Mises criterion can be defined as

(119)

where is the uniaxial yield stress and it is related to the shear yield stress as follows

(120)

and in Eq.(119) is the or defined as

(121)

Concerning the deformation, the elastic-plastic decomposition described in Eq.(115) is rewritten in terms of rates as

(122)

where and are the deformation rates associated to the elastic and plastic responses, respectively.

Combining Eq.(122) with Eq.(98), yields

(123)

The rate of plastic deformations is given by

(124)

where the plastic flow rate is a scalar and represents the plastic flow direction.

Substituting Eq.(124) in Eq.(123), yields

(125)

During plastic loading the plastic flow rate is positive and the state of stress remains on the yield surface . This is consistent with the third Khun-Tucker condition . The consistency condition has the same meaning. Using the chain rule on the consistency condition, yields

(126)

If the yield function depend on the invariant, the following relation holds [IEEEexample:BelyBook, IEEEexample:Prager]

(127)

Combining Eqs.(127) and (125) and substituting in Eq.(126), yields

(128)

In the most plastic models the evolution of the function of the internal variables can be expressed as a function of the plastic strain parameter as follows

(129)

where are the internal variables.

Substituting Eqs.(124) and (129) in 128, the following relation for the plastic strain parameter is obtained

(130)

The plastic flow vector is often derived from a plastic flow potential. If the plastic flow potential coincides with the yield function, the plastic flow is termed . In this case is proportional to the normal of the yield surface, that is . An associative plastic flow has the important advantage that it can lead to a symmetric stiffness matrix [IEEEexample:BelyBook]. In this work an associative plasticity and a constant plastic modulus (for perfect plasticity, ) have been considered. Using these hypotheses the plastic strain parameter is expressed as

(131)

Substituing this relation in Eq.(125), yields

(132)

The same can be computed using a tangent moduli over the whole deformation rate as

(133)

where is the elasto-plastic tangent modulus.

For associative plasticity, the von Mises plastic flow is computed as

(134)

Because the plastic flow vector is deviatoric it follows that

(135)

Form Eqs.(101), (133) and (135), the following expression of the elastoplastic modulus is obtained

(136)

with

(137)

(138)

For perfect plasticity , so and Eq.(136) simplifies to

(139)

Note that the continuum elasto-plastic tangent modulus conserves the symmetry properties of its elastic counterpart. For elastic loading or unloading, .

For a plane strain state, is

In order to guarantee the consistency of the elastoplastic incremental scheme, the so-called algorithm has to be introduced. With this technique, the Khun-Tucker conditions (Eq.(117)) are enforced at the end of a plastic time step in order to recover exactly the yield condition . The return mapping algorithm consists of an initial trial elastic step followed by a plastic corrector one that is activated when the yield function takes a positive value. In Figure 4 from [IEEEexample:Borja], a graphical representation of the return mapping algorithm is shown.
Error creating thumbnail: File missing
Figure 4: Graphical representation of the return mapping algorithm [IEEEexample:Borja].

For associative plasticity, during the plastic corrector step driven by the increment of the plasticity parameter , the plastic flow direction is normal to the yield surface.

For the flow theory and associative plasticity, the return mapping is characterized to be [IEEEexample:SimoTaylor]. This because the von Mises yield surface is circular, thus its normal is also radial. In Figure 5 a graphical representation of the radial return algorithm for plasticity is shown.

Error creating thumbnail: File missing
Figure 5: Graphical representation of the radial return method for plasticity.

The return radial mapping procedure typically starts with the elastic prediction of the stresses. This means that the linear momentum equations are solved using only the elastic part of the continuum tangent modulus (Eq.(101)) and the stress tensor is computed with Eq.(114) (the superindex refers to the iteration index).

Then the effective stress is computed with Eq.(121) and the yield function Eq.(119) is evaluated. If the return mapping iterative procedure is required. In fact, an elastoplastic step has been computed as purely elastic without fulfilling the consistency condition.

The first step of the iterative corrector procedure consists of computing the increment of the plastic parameter as

(140)

For perfect plasticity, the plastic modulus is null, hence the plastic parameter is .

The plastic parameter increment is updated as

(141)

If before this step the material has never suffered plastic deformations, then .

Next the plastic strain and the internal variables are updated according to the plastic correction derived from Eq.(140).

The increment of plastic deformation is

(142)

So the total plastic deformations are

(143)

Once again, for the first plastic step .

Next the deviatoric stresses are updated as

(144)

The plastic flow direction remains unchanged because also the tensor remains unchanged (and radial) during the plastic corrector phase (Eq.(138)). This is a particular feature of the radial return mapping.

From Eqs.(121,144), the updated effective stress is

(145)

Then the yield condition (Eq.119) is verified again with the updated effective stress. If it is not fulfilled, the steps from Eq.(140) to Eq.(144) are repited until .

In Box 5, the return mapping algorithm for the theory and referred to a general elastoplastic time interval is summarized.

2.3.4 Validation examples

In this section several problems are studied in order to validate and the V and the VP elements and to make interesting comparisons. First an example for small displacements is studied. Then a benchmark problem for non-linear solid mechanics, namely the Cook's membrane, is analyzed. The third example is a uniformly loaded circular plate and it involves also plasticity. All these examples are analyzed in statics considering a unique unit time step increment for the velocity-based formulations. The last example is solved in dynamics and for both the hypoelastic and hypoelastic-plastic models.

The first validated problem is a simply supported beam loaded by its self-weight. The problem has been studied in statics so the inertial forces have not been considered. The geometry of the problem is illustrated in Figure 6a and the problem data are given in Table 6.
Error creating thumbnail: File missing
(a) Simply supported beam. Initial geometry.
Figure 6: Simply supported beam. Problem data.

The material properties of the structure can be assimilated to the ones of a structural steel.

The beam undergoes small displacements under the effect of its self-weight, hence linear elastic theory is suitable for computing a reference solution. The accuracy of the formulation is tested by comparing the computed values for the maximum vertical displacement and the maximum XX-component of the Cauchy stress tensor with the values given by a linear elastic analysis. According to this theory, both maximum values are reached in the central section of the beam and they are computed as

(146)

(147)

The problem has been solved in 2D using both the Velocity and the mixed Velocity-Pressure formulations using 3-noded triangular elements. The static problem is solved with only one unit time increment.

In order to verify the convergence of both schemes, the problem has been solved using structured meshes of 3-noded triangles with the following average sizes: 0.25m, 0.125m, 0.05m, 0.025m, 0.0125m. In Figure 7 the finest (mesh size=0.00125m, 32000 elements) and the coarsest (mesh size=0.25m, 80 elements) meshes are illustrated.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
average mesh size= 0.25m
Figure 7: Simply supported beam. Coarsest and finest meshes used for the analysis.
In Figure 8 the solutions for the vertical displacement and the XX-component of the Cauchy stress tensor computed at the Gauss points obtained with the mesh with average size 0.025m are plotted.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Displacements in Y-direction
Figure 8: Simply supported beam. Numerical results.

For the visualization of all the numerical results of this work the pre-postprocessor software GID IEEEexample:GID has been used.

Table 1 collects the values of the maximum vertical displacement (absolute value) and the XX-component of the Cauchy stress tensor computed at the Gauss point using the V-element and the VP-element for different FEM mesh.

0.9

Table. 1 Simply supported beam. Cauchy stress XX-component and maximum vertical displacement for different discretizations.
mesh V-element VP-element
size
2.50E-01 1.46E+06 9.92E-05 1.53E+06 8.41E-05
1.25E-01 2.29E+06 1.37E-04 2.37E+06 1.29E-04
5.00E-02 2.72E+06 1.54E-04 2.80E+06 1.52E-04
2.50E-02 2.82E+06 1.56E-04 2.87E+06 1.56E-04
1.25E-02 2.86E+06 1.57E-04 2.89E+06 1.57E-04

Simply supported beam. Cauchy stress XX-component and maximum vertical displacement for different discretizations. BeamTable

In the examples presented in this section, for the convergence analysis the percentage error is computed versus the solution obtained with the finest discretization as

(148)
For example, the maximum vertical displacement obtained with the finest mesh (average size 0.0125m) is the reference solution for both V and VP elements. The convergence curves for both formulations are plotted in Figure 9. Both elements show a quadratic convergence rate for this error measure.
Simply supported beam. Convergence analysis for the maximum vertical displacement for V and VP elements.
Figure 9: Simply supported beam. Convergence analysis for the maximum vertical displacement for V and VP elements.

The Cook's membrane is a benchmark problem for solid mechanics. The static problem is solved twice in this thesis. In this section a compressible material is considered; in the next chapter the nearly incompressible case is analyzed. In both cases the problem has been solved with only one unit time increment.

The initial geometry of the problem, as well the problem data are given in Figure 10a.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
(a) Initial geometry
Figure 10: Cook's membrane. Initial geometry, material data and FEM mesh (5 subdivisions for each edge corresponding to 50 elements).

The self weight of the membrane has not been taken into account in the analysis, so the membrane deforms under only the effect of the external load applied at its free edge. In this case the structure undergoes large displacements and the solution cannot been computed analytically. The results taken as the reference ones are those published in IEEEexample:Fredriksson. The comparison with the mentioned publication, as well the convergence test are performed for the vertical displacement of point A of Figure 10a with coordinates . According to IEEEexample:Fredriksson, under plane stress conditions this displacement is .

The domain is discretized with a structured mesh and the edges of the membrane have the same number of partitions. In Figure 10 one of the meshes used for this problem is given. The figure refers to the case of 5 elements for each edge of the membrane.

For the 2D simulation a convergence study has been performed using various discretizations. In the finest one the edges have 200 subdivisions, while in the coarsest one only 2.

Error creating thumbnail: File missing
Figure 11: Cook's membrane. Numerical results for the 3D simulation: the XX-component of the Cauchy stress tensor is plotted over the deformed configuration.

The 3D problem (thickness=1) has been solved for an unstructured mesh with average size 0.5 only. The results for this mesh are given in Figure 11, where the XX-component of the Cauchy stress tensor is plotted over the deformed configuration.

For the 3D simulation, the vertical displacements at the central point of the free edge for the V and the VP elements are 23.942 and 23.952, respectively, which correspond to an error versus the reference solution of 0.092% for the V-element and 0.050% for the VP-element.

The vertical displacement of point A of Figure 10a obtained for all the 2D discretizations and for both the Velocity and the Velocity-Pressure formulations is plotted in the graph of Figure 12.
Error creating thumbnail: File missing
Figure 12: Cook's membrane. Vertical displacement of point A of Figure 10a. Results for V and VP elements compared to the reference solution IEEEexample:Fredriksson.

In Table 2 the numerical values are given.

0.9

Table. 2 Cook's membrane. Maximum vertical displacement for different discretizations.
elements V-element VP-element
per edge
2 6.707 7.8105
3 9.0274 10.901
4 11.232 13.515
5 13.1755 15.5985
10 19.037 20.729
15 21.272 22.332
20 22.349 22.987
50 23.658 23.781
100 23.878 23.913
200 23.941 23.95

Cook's membrane. Maximum vertical displacement for different discretizations. CMtable

The convergence curves are given in Figure 13.
Error creating thumbnail: File missing
Figure 13: Cook's membrane. Convergence analysis for the vertical displacement of point A of Figure 10a for the V and VP elements, V-element and VP-element, respectively.

Also in this case the convergence rate is quadratic for both formulations.

The problem analyzed in this section is a simply supported circular plate subjected to a uniform pressure on its top surface. The plate constrains are applied on its lower edge. The plate has a radius and thickness . In this work, the axial symmetry of the problem has not been used, and the plate has been analyzed in 3D using 4-noded tetrahedra. The average size for the tetrahedra is 0.175. This gives 214047 nodes. In Figure 14 the FEM mesh used is shown.
Error creating thumbnail: File missing
Figure 14: Uniformly loaded circular plate. Initial geometry and 3D FEM used.

A hypoelastic-perfectly plastic model has been used, and the problem has been solved with the mixed velocity pressure formulation. For the plastic part, a von Mises yield criterion has been considered. The plate has Young modulus , Poisson ratio and a uniaxial yield stress . The objective of the study is to determine the limit load for the plate. Using the procedure described in IEEEexample:Skrzypek1993, the limit load can be computed analytically by combining the limit analysis and the finite difference method. According to this theory, the limit load can be approximated as

(149)

The same problem was solved in IEEEexample:Peric using eight-noded axisymmetric quadrilateral elements with four Gauss integration points. The limit load obtained with a relatively coarse mesh (10 finite elements distributed in two layers across the thickness) is IEEEexample:Peric.

As in IEEEexample:Peric, the limit load has been considered as the one for which the non-linear procedure cannot longer converge for a small increment of the load.

In Figure 15 the maximum vertical displacement of the plate is plotted against the pressure on the top surface. In Table 3 the numerical values are given.
Error creating thumbnail: File missing
Figure 15: Uniformly loaded circular plate. Maximum deflection versus the applied pressure.

0.9

Table. 3 Uniformly loaded circular plate. Numerical values of the maximum vertical deflection for different applied pressures.
pressure max. deflection pressure max. deflection
101.84 0.0758 260.71 0.677
178.22 0.138 261.21 0.716
229.14 0.236 261.73 0.761
241.87 0.296 262.24 0.816
253.58 0.424 262.73 0.885
258.67 0.567 263.26 0.972
259.69 0.615 263.77 1.088
260.20 0.644 264.27 1.250

Uniformly loaded circular plate. Numerical values of the maximum vertical deflection for different applied pressures. platePericTable

For the present analysis the limit load obtained is , the relative percentage errors with respect the solutions given in IEEEexample:Skrzypek1993 and IEEEexample:Peric are 1.37% and 1.76%, respectively.

In Figure 16 the vertical displacements are depicted over the deformed configuration obtained with the limit load. The plate central section is highlighted in the picture.

In Figures 17 and 18 some snapshots of the von Mises effective stress are plotted over the central section of the plate for the different load conditions. The picture shows clearly the progressive evolution of the plastic zone.
Error creating thumbnail: File missing
Figure 16: Uniformly loaded circular plate. Vertical displacement contours for the maximum pressure sustained by the plate ().
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Overall load=101.84 Overall load=229.14
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 17: Uniformly loaded circular plate. Von Mises effective stress over the deformed configurations for different load conditions (only the central section is depicted) (I/II).
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Overall load=253.58 Overall load=264.27
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 18: Uniformly loaded circular plate. Von Mises effective stress over the deformed configurations for different load conditions (only the central section is depicted) (II/II).

The plane strain cantilever illustrated in Figure 19a has been chosen as the reference case for a large displacement dynamics analysis. The problem data are in given in Table 19. The problem was introduced and studied in IEEEexample:Bely. In the reference publication, the load was applied as a step function at time and its magnitude was defined by the equation where is the coordinate in the direction of the load. However, in this analysis the load has been considered uniformly distributed over the free edge (the overall value is 40, as in IEEEexample:Bely).
Plane strain cantilever. Initial geometry.
(a) Plane strain cantilever. Initial geometry.
Figure 19: Plane strain cantilever. Problem data.

The problem has been solved with both a hypoelastic and a hypoelastic-plastic models. First the results of the hypoelastic model are given.

The problem has been solved in 2D and 3D and using both the V and VP elements. In order to simulate the plane strain state, in the 3D analysis the nodal displacements in the transversal direction to the load have been constrained IEEEexample:Bely.

The reference solution is the elastic one given in IEEEexample:Bely.

For the 2D analysis a convergence study has been performed. Structured finite element meshes have been used and the coarsest and the finest ones have a mean size of 1 and 0.125, respectively. Both meshes are given in Figure 20.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Average mesh size= 1
Figure 20: Plane strain cantilever. Coarsest (mesh size=1, 200 elements) and finest meshes (mesh size=0.125, 12800 triangles) used for the 2D analysis.
For the 3D case, the problem has been solved with the finest mesh only (average size for the 4-noded tetrahedra equal to 0.125). The results for the 3D case obtained with the VP-element are illustrated in Figure 21 where the pressure contours are plotted over the deformed configuration.
Error creating thumbnail: File missing
Figure 21: Plane strain cantilever. Numerical results for the 3D simulation obtained with the VP-element: pressure contours plotted over the deformed configuration.
In Figure 22 the time evolution of the top corner vertical displacement is plotted for each of the FEM meshes. These results have been obtained with the VP-element.
Error creating thumbnail: File missing
Figure 22: Plane strain cantilever. Time evolution of the top corner vertical displacement for different 2D discretizations. Results obtained with the VP-element.

According to IEEEexample:Bely, the maximum vertical displacement is . Table 4 collects the maximum vertical displacement obtained with the V and the VP elements for all the meshes.


Table. 4 Plane strain cantilever. Maximum top corner vertical displacement for different 2D discretizations.
mesh V-element VP-element
size
1 5.759 6.306
0.8 6.144 6.534
0.5 6.568 6.743
0.25 6.811 6.863
0.125 6.875 6.895

Plane strain cantilever. Maximum top corner vertical displacement for different 2D discretizations. BelyTable

The four curves of Figure 23 are the converged time evolution of the top corner vertical displacement obtained with the V-element in 3D, the VP-element in 2D and 3D and the reference solution IEEEexample:Bely. The curves corresponding to the V and VP elements are almost superposed and they match the reference solution.
Error creating thumbnail: File missing
Figure 23: Plane strain cantilever. Time evolution of the top corner vertical displacement. Solutions for the 2D VP-element and the 3D V and VP elements obtained with the finest mesh (average size 0.125) compared to the reference solution IEEEexample:Bely.

The same problem has been solved for an elastic-plastic material with linear hardening. The yield stress is 300 and the plastic modulus is 100. The problem has been solved with the mixed velocity-pressure formulation and by using structured meshes, as the ones of Figure 20. The reference solution is taken from IEEEexample:Bely where the benchmark was proposed. In IEEEexample:Bely the converged value for the maximum top corner vertical displacement is 8.22. The hypoelastic-plastic mixed velocity-pressure formulation converges to 7.97 (error of 2.998% In the graph of Figure 24 the time evolution of the top corner vertical displacement is plotted for the different FEM meshes.
Error creating thumbnail: File missing
Figure 24: Plane strain elastoplastic cantilever. Time evolution of the top corner vertical displacement for different 2D discretizations.

In Table 5 the numerical values for the maximum and the residual top corner vertical displacements are given for each of the FEM mesh.

0.9

Table. 5 Plane strain elastoplastic cantilever. Maximum and residual top corner vertical displacements for different discretizations.
mesh size
1 6.77 3.25
0.8 7.17 3.69
0.5 7.56 4.14
0.25 7.82 4.51
0.125 7.92 4.68
0.1 7.94 4.72
0.0625 7.97 4.77

Plane strain elastoplastic cantilever. Maximum and residual top corner vertical displacements for different discretizations. BelyPlasticTable

The problem has been solved also for the 3D problem for a structured mesh of 4-noded tetrahedra with average size 0.125.

In Figures 25 the von Mises effective stresses are plotted over the deformed configuration at the time instant when the top corner vertical displacement is reached ().
Error creating thumbnail: File missing
Figure 25: Plane strain elastoplastic cantilever. Numerical results for the 3D simulation. Von Mises effective stress plotted over the deformed configuration at ().
In Figure 26 for the same time instant the XX-component of the Cauchy stress tensor is plotted.
Error creating thumbnail: File missing
Figure 26: Plane strain elastoplastic cantilever. Numerical results for the 3D simulation: XX-component of the Cauchy stress tensor plotted over the deformed configuration at ().
In Figure 27 the 3D solution is compared to the 2D results obtained with a structured mesh with the same average size. The curves coincide almost exactly.
Error creating thumbnail: File missing
Figure 27: Plane strain elastoplastic cantilever. Time evolution of the top corner vertical displacement. Numerical results for the 2D and the 3D simulations for the same average mesh size (0.125).

2.4 Summary and conclusions

In this chapter two velocity-based finite element Lagrangian procedures, namely the Velocity and the mixed Velocity-Pressure formulations, have been derived for a general compressible material.

The derivation of both formulations has been carried out with the aim of maintaining the scheme as general as possible. The mixed Velocity-Pressure formulation has been derived exploiting the linearized form of the Velocity formulation. In particular, it has been shown that the tangent matrix of the linear momentum equations is the same for both schemes.

The mixed Velocity-Pressure procedure is based on a two-step Gauss-Seidel solution algorithm. First, the linear momentum equations are solved for the velocity increments, next the continuity equation is solved for the pressure in the updated configuration. At the end of these steps the convergence for the velocities and the pressure is checked. Linear interpolation has been used for both velocity and pressure fields.

Next, both formulations have been particularized for hypoelastic solids. The finite elements generated from the Velocity formulation and the mixed Velocity-Pressure formulations have been called V and VP element, respectively.

The numerical scheme for dealing with associative plasticity has been also given.

Several numerical examples have been given for validating the V and VP elements for large displacements dynamics problems involving both hypoelastic and hypoelastoplastic compressible solids. It has been shown that both elements are convergent for all the numerical examples analyzed.

3 Unified stabilized formulation for quasi-incompressible materials

Chapter 3. Unified stabilized formulation for quasi-incompressible materials

This chapter is devoted to the derivation and validation of the unified stabilized formulation for nearly-incompressible materials. Namely, the cases of quasi-incompressible Newtonian fluids and the quasi-incompressible hypoelastic solids will be analyzed.

Quasi-incompressible materials have a compressibility that is small enough to neglect the variation of density on time but, unlike fully incompressible materials, they are not totally divergence-free and the volumetric strain rate is related to the variation on time of pressure via Eq.(75). This stabilized formulation is based on the mixed velocity-pressure formulation derived in the previous chapter for a general material. In fact a one-field method, as the Velocity formulation presented in Chapter 2, is not sufficient for dealing with the incompressibility constraint. Furthermore the condition IEEEexample:Brezzi imposes the stabilization of the mixed finite element procedure, if an equal order interpolation is used for the velocities and the pressure, as in this work. Consequently, the mixed Velocity-Pressure formulation derived for compressible materials in the previous chapter needs to be stabilized in order to solve quasi-incompressible problems.

In this work a new stabilized Lagrangian method for quasi-incompressible materials is derived. The stabilization procedure is based on the consistent derivation of a residual-based stabilized expression of the mass balance equation using the , also called method IEEEexample:felippaFIC,IEEEexample:FICadt,IEEEexample:Onate,IEEEexample:FICfem,IEEEexample:FICturbolent,IEEEexample:alphaParameter,IEEEexample:FIC2010. The main ideas of this part of the thesis are taken from IEEEexample:nuestroPaper where the stabilization technique was derived for homogeneous viscous fluids. In this chapter it is shown that the procedure can be extended also for the analysis of quasi-incompressible solids.

The FIC approach in mechanics is based on expressing the equations of balance of mass and momentum in a space-time domain of finite size and retaining higher order terms in the Taylor series expansion typically used for expressing the change in the transported variables within the balance domain. In addition to the standard terms of infinitesimal theory, the FIC form of the balance equations contains derivatives of the classical differential equations in mechanics multiplied by characteristic distances in space and time.

In this work the second order FIC form in space and the first order FIC form in time of the mass balance equation have been used as the basis for the derivation of the stabilized formulation. The discretized variational form of the FIC mass balance equation via the FEM introduces terms in the Neumann boundary of the domain and other terms involving the first and second material time derivatives of the pressure. These terms are relevant for ensuring the consistency of the residual formulation.

The FIC stabilization, although is derived using the linear momentum equations, affects only the continuity equation. This means that the general form of the discretized and linearized momentum equations derived in Chapter 2 for the mixed Velocity-Pressure formulation still holds for quasi-incompressible materials. Hence, for hypoelastic quasi-incompressible materials the linear momentum equations are solved through the same linear system derived for the VP (compressible) element in Section 2.3.2. To avoid repetitions, the linearized momentum equations for the mixed velocity-pressure formulation will be recalled only for Newtonian fluids by deriving the tangent matrix and the internal forces according to the constitutive law.

For convenience, the stabilized form of the continuity equation is derived for quasi-incompressible Newtonian fluids, as in IEEEexample:nuestroPaper. Nevertheless it will be shown that the approach can be easily extended to quasi-incompressible hypoelastic solids.

For the fluid analysis a Lagrangian procedure called Particle Finite Element Mehtod (PFEM) IEEEexample:PFEM2006 is used. With the PFEM, the mesh nodes are treated as particles and they move according to the governing equations. The domain is continuously remeshed using a procedure that efficiently combines the Delaunay tessellation and the Alpha Shape Method IEEEexample:Edelsbrunner.

The FIC stabilized formulation here presented IEEEexample:nuestroPaper has excellent mass preservation feature in the analysis of free surface fluid problems. Preservation of mass is a great challenge in the numerical study of flow problems with high values of the bulk modulus that approach the conditions of incompressibility. Mass losses can be induced by the stabilization terms which are typically added to the discretized form of the momentum and mass balance equations in order to account for high convective effects in the Eulerian description of the flow, and to satisfy the condition imposed by the full incompressibility constraint when equal order interpolation of the velocities and the pressure is used in mixed FEMs IEEEexample:BelyBook,IEEEexample:Donea,IEEEexample:Zinc2,IEEEexample:Zinc3.

An important source of mass loss emanates in the numerical solution of free surface flows due, among other reasons, to the inaccuracies in predicting the shape of the free surface during large flow motions IEEEexample:massCons. Mass losses can also occur in the numerical solution of flows with heterogeneous material properties IEEEexample:PFEMmultiF and in homogeneous viscous flows using the Laplace form of the Navier-Stokes equations IEEEexample:NSequa.

In Lagrangian analysis procedures (such as the PFEM) the motion of the fluid particles is tracked during the transient solution. Hence, the convective terms vanish in the momentum equations and no numerical stabilization is needed for treating those terms. Two other sources of mass loss, however, remain in the numerical solution of Lagrangian flows, i.e. that due to the treatment of the incompressibility constraint by a stabilized numerical method, and that induced by the inaccuracies in tracking the flow particles and, in particular, the free surface.

The discretized variational form of the FIC mass balance equation via the FEM introduces terms in the Neumann boundary of the domain, and other terms involving the first and second material time derivatives of the pressure that are relevant for ensuring the consistency of the residual formulation. These terms are also crucial for preserving the mass during the transient solution of free surface Lagrangian flows. In addition they enable the computation of the nodal pressures from the stabilized mass balance equation without imposing any condition on the pressure at the free surface nodes, thus eliminating another source of mass loss which occurs when the pressure is prescribed to a zero value on the free surface in viscous flows.

A section of this chapter is exclusively dedicated to show the excellent mass preservation features of the PFEM-FIC stabilized formulation for free surface flow problems.

Various approaches have been developed in the recent years for approximating fluid flows by means of a quasi-incompressible material. In practice, they all consider the Navier-Stokes problem with a modified mass conservation equation where a slight compressibility is added to the fluid. In previous works, see for example IEEEexample:alphaParameter,IEEEexample:Turkel,IEEEexample:Chorin,IEEEexample:Zinc3, the fluid compressibility was introduced by relaxing the incompressibility constraint by means of a penalty parameter. Alternatively, the small compressibility of the fluid can be introduced by considering the actual bulk modulus of the fluid , which gives this operation a physical meaning, IEEEexample:Donea,IEEEexample:Idelshon,IEEEexample:Zinc3. The success of quasi-incompressible formulations in fluid mechanics relies on their important advantages from the numerical point of view. The most obvious one is that the quasi-incompressible form of the continuity equation yields a direct relation between the two unknown fields of the Navier-Stokes problem, the velocities and the pressure (Eq.(75)). This is useful if the problem is solved using a partitioned scheme because the velocity-pressure relation is crucial for deriving the tangent matrix of the momentum equations. Furthermore, another important drawback of fully incompressible schemes is eluded; the incompressibility constraint leads to a diagonal block of a zero matrix in the global matrix system. Consequently, a pivoting procedure is required to solve numerically this kind of linear system. It is well known that the computational cost associated to this operation is high and it increases with the number of degrees of freedom of the problem. The compressibility terms that emanate from the quasi-incompressible form of the continuity equation fill the diagonal of the global matrix, thus overcoming these numerical difficulties.

On the other hand, quasi-incompressible schemes insert in the numerical model parameters that have typically high values and can lead to different numerical instabilities. For example, large values of the penalty parameter or, equally, physical values of the bulk modulus, can compromise the quality of the analyses, or even prevent the convergence of the solution scheme, IEEEexample:Zinc3. For this reason, generally, the value of the actual bulk modulus is reduced arbitrarily to the so-called . Nevertheless, an excessively small value of the pseudo bulk modulus changes drastically the meaning of the continuity equation of the original Navier-Stokes problem. In other words, the incompressibility constraint would not be satisfied at all. Furthermore, the bulk modulus is proportional to the speed of sound propagating through the material. Hence, we have to guarantee that the order of magnitude of the velocities of the problem is several times smaller than the velocity of sound in the medium.

In this work, the pseudo-bulk modulus is used for the tangent matrix of the linear momentum equations, while the actual physical value of the bulk modulus is used for the numerical solution of the mass conservation equation. The pseudo bulk modulus is computed as a proportion of the real bulk modulus of the fluid through the parameter such that with . A new numerical strategy for computing a priori the optimum value for the pseudo bulk modulus is also derived. For free surface flow problems, it will be shown that the scheme guarantees the good conditioning of the linear system, good convergence and excellent mass preservation features.

The lay-out of this chapter is the following. First the stabilized FIC form of the mass balance equation is derived. The procedure is described for a Newtonian fluid from the local and continuous form until the fully discretized matrix form. Next the linearized and discretized expression of the linear momentum equations derived in Chapter 2 for the mixed Velocity-Pressure formulation is adapted for Newtonian fluids and the complete solution scheme is given. Next the FIC stabilization procedure is extended to hypoelastic quasi-incompressible materials and the complete solution scheme is given also for this constitutive model.

In Section 3.4, free surface flows are analyzed in detail. First the essential features of the PFEM are given, then the mass preservation feature of the PFEM-FIC stabilized formulation is shown together with some explicative numerical examples. Finally the conditioning of the linear system is studied and a practical and efficient technique to predict the optimum value for the pseudo bulk modulus is given.

The chapter ends up with several validation examples for quasi-incompressible solid and fluid mechanics problems.

3.1 Stabilized FIC form of the mass balance equation

The FIC stabilized form of the continuity equation is here derived. For convenience, the derivation procedure is carried out for quasi-incompressible Newtonian fluids.

3.1.1 Governing equations

For the sake of clarity, the local forms of the linear momentum and the continuity equations are recalled. From Eqs.(1) and (75) yields

(150)

(151)

The terms of Eqs.(150 and 151) have already been defined in the previous chapters.

The standard form of the constitutive relation for a Newtonian fluid reads

(152)

where is the viscosity and the deviatoric strain rate is defined from Eq.(71) as

(153)

where is the volumetric strain rate.

Substituting Eqs.(152) and (153) into (150), gives a useful form of the momentum equations for the Newtonian fluids as

(154)

3.1.2 FIC mass balance equation in space and in time

Previous stabilized FEM formulations for quasi and fully incompressible fluids and solids were based on the first order form of the Finite Calculus (FIC) balance equation in space IEEEexample:FICadt–IEEEexample:FICfsi2,IEEEexample:FIC2,IEEEexample:FICfem–IEEEexample:FICturbolent,IEEEexample:alphaParameter. In this work, for the derivation of stabilized formulation both the second order FIC form of the mass balance equation in space IEEEexample:alphaParameter,IEEEexample:FIC2010 and the first order FIC form of the mass balance equation in time are used. These forms read respectively IEEEexample:nuestroPaper

(155)

and

(156)

Eq.(155) is obtained by expressing the balance of mass in a rectangular domain of finite size with dimensions (for 2D problems), where are arbitrary distances, and retaining up to third order terms in the Taylor series expansions used for expressing the change of mass within the balance domain. The derivation of Eq.(155) for 2D incompressible flows can be found in IEEEexample:FIC2010.

Eq.(156), on the other hand, is obtained by expressing the balance of mass in a space-time domain of infinitesimal length in space and finite dimension in time IEEEexample:FICadt.

The FIC terms in Eqs.(155) and (156) play the role of space and time stabilization terms respectively. In the discretized problem, the space dimensions and the time dimension are related to characteristic element dimensions and the time step increment, respectively as it will be explained later.

Note that for and the standard form of the mass balance equation (151), as given by the infinitesimal theory, is recovered.

3.1.3 FIC stabilized local form of the mass balance equation

Substituting Eq.(151) into Eqs.(155) and (156) the second order FIC form in space and the first order FIC form in time of the mass balance equation for a general quasi-incompressible material read

(157)

and

(158)

The FIC form of the mass balance equation is expressed in terms of the momentum equations. Neglecting the space changes of the viscosity , from Eq.(154) the following expression is obtained

(159)

Hence

(160)

In the above two equations is a static momentum term defined as

(161)

Substituting Eq.(160) into Eq.(157) and neglecting the space changes of and in the derivatives, the following form is obtained

(162)

Observation of the term involving the material derivative of in Eq.(162) gives

(163)

Substituting Eq.(163) into (162) gives

(164)

From Eq.(158),

(165)

Substituting Eq.(165) into Eq.(164) gives

(166)

Multiplying Eq.(166) by gives, after grouping some terms,

(167)

After some further transformations,

(168)

Where is a stabilization parameter given by

(169)

For transient problems the stabilization parameter is computed as

(170)

where is the time step used for the transient solution and is a characteristic element length.

The coefficient multiplying the second space derivatives of in Eq.(168) is much smaller than the coefficients multiplying the rest of the terms in this equation. Numerical tests have shown that the results are not affected by this term. Consequently, this second space derivative term will be neglected in the rest of this work.

Hence the FIC stabilized form of the mass balance equation is written as

(171)

Note that the term within in Eq.(171) (see the definition of in Eq.(161)) vanishes for a linear approximation of the velocity field. This is the case for the simplicial elements used in this work.

3.1.4 Variational form

Multiplying Eq.(171) by arbitrary (continuous) test functions (with dimensions of pressure) and integrating over the analysis domain gives

(172)

Integrating by parts the last integral in Eq.(172) (and neglecting the space changes of ) yields

(173)

where are the components of the unit normal vector to the external boundary of .

Using Eq.(160) an equivalent form for the boundary term BT of Eq.(173) is obtained as

(174)

where is the derivative of the volumetric strain rate in the direction of the normal to the external boundary and is the velocity normal to the boundary.

The term can be approximated as follows

(175)

where and are respectively the values of and at exterior and interior points of the boundary and is a characteristic length in the normal direction to the boundary. Figure 28 shows an example of the computation of at the side of a 3-noded triangle adjacent to the external boundary. The same procedure applies for 4-noded tetrahedra.

Computation of the term  of \displaystyle \frac23 μퟃdv\over ퟃn at the side ij of a 3-noded triangle ijk adjacent to the external boundary Γ.
Figure 28: Computation of the term of at the side of a 3-noded triangle adjacent to the external boundary .

Clearly, at external boundaries and . Hence, coincides with the volumetric strain in the 3-noded triangular element adjacent to the boundary.

Using above argument Eq.(175) simplifies to

(176)

On the other hand, the stresses at any boundary satisfy the traction equilibrium condition

(177)

Substituting Eqs.(152) and (153) into (177) and multiplying all terms by , yields

(178)

where is the normal traction to the boundary (Figure 28) and .

From Eq.(178)

(179)

Substituting Eq.(179) into Eq.(176) and this one into Eq.(174), gives the expression of the boundary integral of Eq.(173) as

(180)

The normal velocity is fixed at a Dirichlet boundary and, hence, at . Also, accepting that at , the surface tractions at coincide precisely with the reactions computed as . Hence, the boundary integral can be neglected at a Dirichlet boundary and, therefore, it has a meaning at a Neumann boundary only. In conclusion,

(181)

Substituting Eq.(181) into (173) and using the expression of of Eq.(161) yields the variational expression of the stabilized mass balance equation, after rearranging the different terms, as

(182)

Expression (182) holds for 2D and 3D problems.

The terms involving the first and second material time derivative of the pressure and the boundary term in Eq.(182) are important to preserve the consistency of the residual form of the FIC mass balance equation. This, in turn, is essential for preserving the conservation of mass in the transient solution of free flow problems. The form of Eq.(182) is a key contribution of the new FIC-based stabilized formulation, versus previous works on this topic IEEEexample:FIC1,IEEEexample:FICale,IEEEexample:alphaParameter,IEEEexample:FIC2010,IEEEexample:PavelPaper.

At an unloaded free surface (Neumann) boundary , and hence

(183)

For an inviscid fluid and Eq.(183) simplifies to

(184)

Accounting for the term in the boundary integral of Eqs.(182-184) has proven to be relevant for the enhanced conservation of mass in free surface flows IEEEexample:nuestroPaper. On the other hand, the effect of the term involving was negligible in all the problems solved in this work IEEEexample:nuestroPaper.

Eq.(182) is the starting point for deriving a new class of linear triangles with discontinuous pressure field adequate for analysis of incompressible flows with heterogeneous material properties IEEEexample:onateP1P2.

The discretization of the mass stabilized equation is performed as it has been shown for the continuity equation in the previous chapter. So the analysis domain into finite elements is discretized using 3-noded linear triangles () for 2D problems and 4-noded tetrahedra () for 3D problems with local linear shape functions defined for each node () of an element .

3.1.5 FEM discretization and matrix form

Substituting the approximations (9) and (77) into Eq.(182) and choosing a Galerkin form with gives the discretized form of the stabilized mass balance equation, after eliminating the arbitrary test functions as

(185)

where

(196)

Eqs.(185) can be written in matrix form as

(197)

The matrices and vectors in Eqs.(197) are assembled from the element contributions as

(198)

(199)

(200)

(201)

(202)

(203)

The boundary terms in vector (Eq.(203)) can be incorporated with the matrices of Eq.(197). This, however, leads to a not symmetrical set of equations. For this reason these boundary terms are computed iteratively within the incremental solution scheme.

In a free surface fluid the presence in Eq.(197) of matrix (Eq.(200)) enables the computation of the pressure without the need of prescribing its value at the free surface. This eliminates the error introduced when the pressure is prescribed to zero in free boundaries, which leads to considerable mass losses for viscous flows. Matrix was introduced into the discretized stabilized mass balance equation in IEEEexample:massCons using a fractional step method and heuristic arguments.

3.2 Solution scheme for quasi-incompressible Newtonian fluids

The solution scheme for quasi-incompressible Newtonian fluids has a structure very similar to the two-step procedure presented for a general compressible material in Chapter 2. The linear momentum equations are solved for the increment of the velocities and the pressures are obtained from the continuity equation in the updated configuration. However, for dealing with the incompressibility, the stabilized form of the continuity equation (Eq.(197)) has to be considered. Furthermore, the linearized form of the momentum equations (Eq.(64)) has to be modified according to the constitutive law of Newtonian fluids.

3.2.1 Governing equations

For the sake of clarity, the general linearized form of the momentum equations is recalled. For each iteration the following linear system is solved

(204)

where:

(205)

with

(206)

(207)

(208)

and

(209)

In order to obtain the linearized form of the momentum equations for quasi-incompressible Newtonian fluids, the tangent constitutive tensor in matrix (206) and the Cauchy stress tensor that appears in the residual vector and in have to be computed using the adequate constitutive law.

For a Newtonian fluid, the stresses can be computed as

(210)

For quasi-incompressible materials Eq.(151) holds within a general time interval . For clarity purposes, the relation is rewritten as

(211)

Considering a time interval and substituting Eq.(211) into (210) yields

(212)

where the fourth-order tensor has been defined in Eq.(97).

For convenience, Eq.(212) is rewritten in the following form

(213)

where the following substitutions have been done:

(214)
(215)
(216)

The goal is to obtain a relationship between the measure of rate of stress and the rate of deformation in the form of Eq.(43). Eq.(213) shows that, according to the constitutive law for Newtonian fluids, the rate of deformation is related to the Cauchy stress and not to rate of the Cauchy stress, as for hypoelastic solids (Eq.(98)). For this reason, in fluids preserving the objectivity of the stress rate measures is not a critical issue as for hypoelastic solids. Rigid rotations do not cause any state of stress, because the Cauchy stress is directly obtained from the rate of deformation. For these reasons, the rate of Cauchy stress can be simply defined from Eq.(213) as

(217)

Note that Eq.(217) has the same structure as Eq.(43). For convenience, the Newtonian tangent moduli for the rate of stress is computed as the sum of the deviatoric and volumetric parts. Hence, substituting into Eq.(206) yields that for a quasi-incompressible Newtonian fluid the material part of the tangent matrix can be written as

(218)

where and are defined for a generic finite element and the pair of nodes as

(219)

(220)

The volumetric part of the tangent matrix can compromise the conditioning of the linear system because its terms are orders of magnitude larger than the viscous part. In order to prevent the numerical instabilities originated by the ill-conditioning of the tangent matrix, a reduced pseudo-bulk modulus, computed from the actual bulk modulus as , can be used in the expression of without altering the numerical results IEEEexample:nuestroPaperBulk. An adequate selection of the pseudo-bulk modulus also improves the overall accuracy of the numerical solution and the preservation of mass for large time steps IEEEexample:nuestroPaperBulk. For fully incompressible fluids (), a finite value of is used in as this helps to obtaining an accurate solution for velocities and pressure with reduced mass loss in few iterations per time step IEEEexample:nuestroPaperBulk. These considerations, however, do not affect the value of within matrix in Eq.(197). Clearly, the value of the terms of can also be limited by reducing the time step size. This, however, increases the overall computational cost. Another approach for improving mass conservation in incompressible flows was proposed in IEEEexample:PavelPaper. In Section 3.4.3 the details about this procedure will be given.

The pressure is obtained from the stabilized form of the mass balance equation (Eq.(197)). Introducing the time integration of the pressure (Eqs.(79, 80)) into Eq.(197), yields

(221)

where

(222)

and

(223)

where all the matrices and the vectors of Eqs.(221-223) have been defined in Eqs.(198-203)

3.2.2 Solution scheme

The complete solution scheme for a quasi-incompressible Newtonian fluid is described for a generic time interval in Box 8.

In Box 9 all matrices and vectors that appear in Box 8 are given.

3.3 Solution scheme for quasi-incompressible hypoelastic solids

The differences between the solution scheme for quasi-incompressible and compressible hypoelastic solids concern only the continuity equation. In fact, the same linearized form of the continuum equations derived in the previous chapter for the compressible hypoelastic model within the mixed Velocity-Pressure formulation holds also for the quasi-incompressible scheme. However, for dealing with incompressibility the scheme needs to be stabilized

In this work the FIC stabilization procedure originally derived for Newtonian fluids in IEEEexample:nuestroPaper and proposed again in Section 3.1 of this chapter, has been tested for the first time for quasi-incompressible hypoelastic solids. In order to use the same form of Eqs.(221-223) for hypoelastic quasi-incompressible solids, the fluid parameters (the fluid viscosity and the fluid bulk modulus ) should be substituted with the equivalent solid parameters. The similarity between Newtonian fluids and hypoelastic solids is evident comparing the computation of the Cauchy stress tensor increment for the two mentioned constitutive laws.

For quasi-incompressible Newtonian fluids Eq.(213) holds and, for clarity purposes, is here rewritten as

(224)

From Eqs.(111), the increment of the Cauchy stress for hypoelastic solids is

(225)

where bulk modulus for the solid is computed from the Lamè constants and or in terms of the Young modulus and the Poisson ratio of the material (Eqs.(92, 93,96)).

Eqs.(224 and 225) show the duality between hypoelastic and Newtonian quasi-incompressible constitutive laws. In the latter the deviatoric and the volumetric parts of the Cauchy stress tensor are controlled by the dynamic viscosity and the bulk modulus , respectively. The equivalent roles in hypoelastic solids are taken by the second Lamè constant scaled by the time increment () and the bulk modulus . Thanks to this equivalence, the stabilized mass continuity equation derived for quasi-incompressible fluids (Eq.(197)) holds also for quasi-incompressible hypoelastic solids just by replacing the fluid parameters and with the equivalent terms for hypoelastic solids and .

Note that using this analogy the differences between the incremental solution schemes for Newtonian fluids and quasi-incompressible hypoelastic solids are minimal and they essentially differ on the way to compute the stresses and the tangent moduli according to the specific constitutive law. For hypoelastic quasi-incompressible model, the stress tensor is computed with Eq.(114) and the tangent moduli is obtained from Eq.(101). The solid finite element implemented with this formulation is called VPS-element.

In Box 10, the iterative solution incremental scheme for quasi-incompressible hypoelastic solids using the stabilized mixed velocity-pressure formulation is described for a generic time increment .

In Box 11 all matrices and vectors that appear in Box 10 are given.

3.4 Free surface flow analysis

Free surface flows are those fluids with at least one side in contact with the air. Their numerical solution is critical because the free surface contours change continuously and they need to be tracked in order to solve accurately the boundary value problem. In this work, this task is carried out using the Lagrangian finite element procedure called IEEEexample:CIMNEPFEM. In the PFEM the mesh nodes are treated as the fluid particles. As a consequence, the free surface contours are automatically detected by the nodes positions. The PFEM will be described and analyzed in detail in Section 3.4.1.

Another typical drawback associated to free surface flows is the deterioration of the mass preservation of the numerical method. The indetermination of the free surface position introduces in the scheme an additional mass loss source to those induced by the inaccuracy of the numerical method. In Section 3.4.2 the mass conservation feature of the PFEM-FIC stabilized formulation is studied in detail and many validation examples are given.

In Section 3.4.3 the key role of the pseudo bulk modulus for guaranteeing the mass preservation of the quasi-incompressible formulation will be highlighted. In fact, the actual value of the bulk modulus deteriorates the conditioning of the linear system. This may affect the mass preservation feature and even the convergence of the whole scheme. Using a reduced value for the bulk modulus, computed as , all these inconveniences are overcome. In the same section a practical strategy for computing the pseudo bulk modulus is also given and tested with several numerical examples.

3.4.1 The Particle Finite Element Method

The is a Lagrangian finite element procedure ideated and developed by Idelsohn and Oñate and their work group IEEEexample:PFEM2004,IEEEexample:PFEM2006,IEEEexample:PFEM2008,IEEEexample:CIMNEPFEM. The PFEM addresses those problems where severe changes of topology occur. This may happen in free surface fluid dynamics, in non-linear solid mechanics, in FSI problems or in thermal coupled problems. The essential idea of the PFEM is to follow the topology of the deformed bodies by moving the nodes of the mesh according to the equations of motion in a Lagrangian way. In other words, in the PFEM the mesh nodes are treated as particles and they transport their momentum together with all their physical properties. All this produces the deformation of the finite elements discretization that needs to be rebuilt whenever a threshold value for the distortion is reached. Remeshing is one of the most characteristic points of the PFEM IEEEexample:PFEMmesh. This operation is performed via an efficient combination of the Delaunay Tessellation IEEEexample:SaalfeldDelaunay,IEEEexample:EdelsbrunnerDelaunay and the Alpha Shape method IEEEexample:Edelsbrunner. Once the mesh is generated, the differential problem is integrated again over the new mesh in the classical FEM fashion.

In the following parts of this section the essential points of the PFEM will be highlighted. Specifically, first the remeshing procedure will be described, then the basic steps of the PFEM will be summarized and finally the advantages and disadvantages of the technique will be highlighted.

3.4.1.1 Remeshing

In order to solve accurately the FEM problem, at each time step it has to be guaranteed a good quality for the mesh. In fact, even a single degenerated element may compromise the entire computation. For this reason the remeshing algorithm must be reliable and robust. In the PFEM this is guaranteed through an efficient algorithm that combines two trustworthy techniques, the Delaunay triangulation and the Alpha Shape method.

Given a cloud of points in the space, the Delaunay tessellation guarantees the creation of the most homogenous discretization for those points. This is obtained by ensuring that the circumball, cirmumcircle in 2D, built on all the nodes of a simplex does not contain any other node of some other simplex IEEEexample:CohenDelaunay. Actually this is guaranteed only in 2D where it is proved that there exists a unique Delaunay triangulation. In 3D there may be the eventuality that the criterion is not fulfilled for certain tetrahedra. In 2D the Delaunay triangulation has the property of maximizing the minimum inner angle of all the triangles of the mesh (this is not guaranteed in 3D). Note that the maximization of the minimum angle does not imply the minimization of the maximum angle, see Figure 29.
Set of pointsDelaunay Not  Delaunay
Set of points Not Delaunay
2D Delaunay triangulation of a set of points.
Figure 29: 2D Delaunay triangulation of a set of points.

This strategy may be helpful for avoiding the creation of degenerated elements as the (simplices with null surface or volume) but it is not sufficient for their complete elimination. It will be explained later that further controls and modifications on the cloud of points are required for this objective.

The Delaunay triangulation alone cannot ensure the detection of the physical boundaries of the bodies. For this purpose, the so-called Alpha Shape method is applied on the tessellation obtained with the Delaunay procedure. The role of the Alpha Shape method is to erase those simplices that are excessively distorted or big. In order to do this, the mean mesh size and the circumradius of each element are computed. Next the following check is done for all the elements

(226)

where is a parameter that is chosen arbitraly.

The selection of is arbitrary. However, it must be taken into account the following considerations. If is excessively large, the Alpha Shape method is worthless because no element will be erased. This means that the Delaunay triangulation remains unchanged. On the other hand, an excessive small value for induces the elimination of all the simplices. Specifically the lowest admissible parameter is the one of the equilateral triangle for which . In Figure 30 the circumcircle is plotted for simplices with the same mean length for edges. Given a mean size for the edges of a simplex, the equilateral triangle is the one with the minimum circumradius. Equally, for a given circumradius, the equilateral triangle is the simplex which edges have the maximum mean length.
Equilateral triangleSharp triangle Flat triangle
(a) Equilateral triangle (b) Flat triangle
Error creating thumbnail: File missing
Figure 30: Triangles with the same mean size and their circumcircles.
The arbitrariness in the choice of the parameter may affect the numerical results. In fact, two different values for may give two different configurations. See for example the case of Figure 31. For the same Delaunay triangulation (Figure 31b), the configuration obtained with (Figure 31c) is different from the one obtained with (Figure 31d), where .
Set of pointsDelaunay triangulation Alpha Shape (α₁)
(a) Set of points (b) Alpha Shape ()
Alpha Shape (α₂) 2D triangulation of a set of points. Alpha shape method applied on the Delaunay triangulation (α₁> α₂).
(c) Alpha Shape ()
Figure 31: 2D triangulation of a set of points. Alpha shape method applied on the Delaunay triangulation ().

Specifically, the mesh obtained with accepts a larger number of elements from the Delaunay discretization than the one given by , especially in the free surface zone. However, in the section dedicated to the validation examples it will be shown that the method is convergent. This means that the error introduced by a certain reduces with the refinement of the mesh.

The remeshing can be performed also in the opposite order, hence using first an Alpha Shape method for detecting the contours of the domains and then performing the Delaunay triangulation with the geometric restrictions of those boundaries IEEEexample:PFEM2008,IEEEexample:PFEMzhang1,IEEEexample:PFEMbederosion. This can be done by performing the Delaunay triangulation via the technique IEEEexample:AdvancingFront.

Both remeshing algorithms may generate particles or elements separated from the rest of the domain. Remember that all the information is stored in the nodes, so the governing equations can be computed also for an isolated particle. This allows, for example, the simulation of detached fluid drops.

However, the first strategy is preferable because it leads to a reduced computational cost and it requires a lower implementation effort. For these reasons, in this work the remeshing is performed applying the Alpha Shape technique after the Delaunay triangulation.

Neither the Alpha Shape method nor the Delaunay triangulation can prevent from the creation of highly distorted elements. Let consider, for example, the triangle depicted in Figure 30b. The circumradius associated to this type of triangle is not excessively large and it may fulfill Eq.(226). The Alpha Shape method and the Delaunay triangulation can only ensure the best physical mesh (so the best discretization that respects the physical boundaries of the domain) for a given cloud of points. However, during the motion, the nodes distribution may be such that it is impossible to avoid the presence of poor quality elements. In order to avoid this eventuality the remeshing with the Delaunay tessellation and the Alpha Shape method must be combined with some additional controls.

For example, for avoiding the presence of sharp elements in the triangulation, as the one shown in Figure 32, one should monitor the edge lengths of each triangle.
Distorted element in a 2D triangulation.
Figure 32: Distorted element in a 2D triangulation.

Whenever the distance between two nodes is smaller than a prearranged critical length, one of those nodes is removed and placed in other zone of the mesh, for example where there are some distorted elements. The reallocation of the nodes is done with the purpose of preserving the number of particles and, consequently, guaranteeing the conservations of the mean mesh size (if the volume does not change).

Another type of critical elements are the slivers, so the simplices with null area or volume. These may form when there are three points contained in the same line, or four points laying in the same 3D plane as shown in Figure 33.
Sliver element in 3D.
Figure 33: Sliver element in 3D.

In theory, the Alpha Shape method should not allow the formation of those elements because the associated circumradius is huge, see Figure 30c. However, these elements are critical from the numerical point of view and they can escape to the Alpha Shape control. Furthermore, also in the case that the sliver element is erased from the Delaunay triangulation, this would create a not physical void within the domain. For all these reasons, the formation of sliver elements must be prevented. For these elements the control for the minimum edges lengths is worthless because the sliver may form also when the element nodes are not close. In this work, the formation of slivers is prevented by checking the element areas. Whenever the element surface or volume is less than a threshold value, computed for example with respect to the mean area of the mesh, one of the nodes of those elements is removed. Note that this control is also beneficious versus the formation of sharp elements.

After these operations over the nodes, the discretization needs to be updated. This can be done modifying locally the mesh or performing again the Delaunay triangulation and the Alpha Shape control. Once the updated mesh is created, the nodal variables of the new particles are assigned, via interpolation, with the values of the neighbor nodes.

All these criteria have been explained for homogeneous meshes. However, they can be easily extended to heterogeneous meshes. For example, if there is a mesh refinement in a specific zone of the domain, that refined mesh can be preserved if the described checks for the minimum nodes distance and the elements surface are done not for the whole mesh but for the local mesh.

3.4.1.2 Basic steps

The basic steps of the PFEM algorithm are here given with the help of a graphic representation (Figure 34).
Error creating thumbnail: File missing
Figure 34: Sequence of steps to update a “cloud” of fluid nodes and a discretized solid domain from time () to time ().

Consider a domain containing fluid and solid subdomains.

In this work the PFEM is used only for the fluid domains. The solids are computed with the classical FEM and they maintain the same discretization for the whole duration of the analysis. In many other works the PFEM has been also used for modeling the solid and the fluid domains, indifferently IEEEexample:PFEMexcavation,IEEEexample:PFEMtunnelling,IEEEexample:PFEMmelting,IEEEexample:OliverPFEMsolid.

At the beginning of a time step (first picture of Figure 34) the fluid is represented by a cloud of points, or particles, that store all the information about the physics (for example the density and the viscosity of the fluid), the geometry (information about the local mesh size), the kinematics (the velocity and the acceleration) and the pressure. On the other hand, the solid has the same mesh of the previous time step and the history for its kinematics and the stress and strain fields is stored in both the nodes and the elements.

The first step represents the distinguishing step of the PFEM and it consists on the creation of a mesh for the fluid. The algorithm is the one that has been described in Section 3.4.1.1. So first the Delaunay triangulation is performed and then the physical contours of the domains are recognized with the Alpha Shape method. Specifically, the fluid detects automatically all its boundaries, hence the rigid walls, its free surface contours and the interface with the solid. During this step flying subdomains may form or/and some particles may detach from the rest of the domain (second picture of Figure 34)

Once all the domains are discretized, it is possible to solve the equations of continuum mechanics for each of the subdomains using the FEM. The state variables are computed at the next (updated) configuration for : velocities, pressure and viscous stresses in the fluid and displacements, stresses and strains in the solid. After that, the fluid mesh can be erased and these steps are repeated. For reducing the computational cost associated to the detection of the fluid boundaries it is very useful to memorize those nodes that in the previous time step were located at the contours. In this way, the nodes checked by the Alpha Shape method can be drastically reduced.

Once again, note that in the fluid all the information is stored in the nodes so this operation does not require interpolation procedures for recovering the elemental information. If the PFEM would be used also for the solid mechanics analysis, an algorithm for transferring the elemental information from the previous to the new mesh should be implemented. In IEEEexample:CarbonellThesis an efficient technique for this purpose is described and validated.

3.4.1.3 Advantages and disadvantages

The key differences between the PFEM and the classical FEM are the remeshing technique and the identification of the domain boundary at each time step. The quality of the numerical solution depends on the discretization chosen as in the standard FEM and, as it has been explained, adaptive mesh refinement techniques can be used to improve the solution. So from this point of view, during the computation step, the PFEM behaves as a classical Lagrangian FEM, with the same advantages and drawbacks.

However with the PFEM the remeshing may introduce a further source of error in the numerical scheme. In fact during this operation some new elements may be created and some other erased. These changes cause local perturbations of the equilibrium reached at the previous time step and they may affect the accuracy and the convergence of the scheme. In fact at the beginning of the time step, each node stores the kinematics and the pressure obtained with a discretization that may be different than the new one. So those values that at the end of the previous step were at equilibrium for the previous configuration and mesh, may not be at equilibrium at the beginning of the new time step for the new configuration created by the remeshing.

Another drawback induced by the remeshing concerns the variation of volume. Although, generally, the volumes gained and lost with the remeshing are compensated, there may be cases where the variation of volume is not negligible. From the personal experience of the author, remeshing tends to increase the volume of the domain. For avoiding this inconvenience, some additional criteria for limiting the creation of new elements should be added. Generally, the critical zones are the ones near to the free surface, where, typically, the remeshing tends to create new elements that 'close' the surface waves and increase the overall volume of the fluid, as shown in Figure 35.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
(a) Mesh at the end of step (b) Mesh detail at the end of

step

Error creating thumbnail: File missing
Error creating thumbnail: File missing
(c) Mesh detail at the beginning of step +1
Figure 35: Example of variation of volume induced by the remeshing.

It is possible to reduce this tendency of the remeshing by locally penalizing the Alpha Shape method. In this work, the free surface nodes have assigned a smaller Alpha Shape parameter than the other nodes of the mesh.

This penalization is also good for improving the timing of the contact with the rigid walls or the solid interfaces. In particular, this strategy may delay the phenomenon illustrated in Figure 36. The elements created at the free surface may accelerate the impact between the fluid streams and the containing walls.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
End of time step
Figure 36: Example of detection of contact through the remeshing.

The smaller is for the interface nodes, the later a contact element between those nodes and the solid or rigid contour is built. In fact, a smaller for the interface nodes delays the formation of the contact elements with the rigid (or deformable) contours reproducing better the physical phenomenon of the contact. However, the Alpha Shape cannot be penalized excessively because otherwise the fluid could pass through the wall, or the solid interface.

The remeshing also represents an additional computational cost. In previous works it has been found that the computational time associated to the remeshing grows linearly with the number of nodes IEEEexample:PFEMfsi. Specifically, for a single processor Pentium IV PC the meshing consumes for 3D problems around 15% of the total CPU time per time step, while the solution of the equations (with typically 3 iterations per time step) and the system assembly consumes approximately 70% and 15% of the CPU time per time step, respectively. The treatment of the boundary nodes is an issue that deserves a particular attention. In the previous works with the PFEM stick conditions were generally used for the nodes on the rigid walls. However, this may induce pressure concentrations on the boundary elements and deteriorate locally the quality of the mesh. In fact the stick condition affects a zone that has the same order of magnitude of the spatial discretization. This may induce a non-physical behavior, especially when inviscid fluids are analyzed or/and coarse meshes are employed. The problem is even more critical in 3D. Recently, Cremonesi IEEEexample:CremonesiSlip2,IEEEexample:CremonesiSlip1 used slip conditions in the simulation of landslides for better modeling the interaction between a landslide and the substrate interface. In the mentioned works, the wall nodes are treated in an Eulerian way by adding the convective term for the boundary elements. In this work, the wall particles are still computed in a Lagrangian framework and they are free to move along the direction of the walls. The slip conditions are simulated with a simple algorithm. Essentially it consists on leaving the wall particles move along the direction of the wall until when the separation from the original position is larger than a prearranged critical distance. In that case, during the remeshing procedure, the particle is removed and reallocated at its original position. The reallocation of the particle is done in order to prevent the creation of voids in the walls that may cause fluid leakage. The kinematics and the pressure of the moved particle are obtained via interpolation from the neighbor nodes, in the classical PFEM fashion. In Figure 37 an example of the application of the algorithm is given.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 37: Agorithm for moving the wall particles.

The pictures refer to the front of a flow.

The void circles follow the position of two nodes of the wall discretization, the full circles denote the original position of those nodes. The pictures show that the wall nodes are free to move along the wall direction and when they reach the maximum separation from the original position (in this case , where is the distance between the node of the wall) they are reallocated at the original position. In the section dedicated to the validation examples, some interesting comparisons between stick and slip problems are given, both in 2D and 3D. It will be shown that this algorithm can be very helpful for preserving the quality of the mesh and for the overall accuracy of the computation.

All the weak points of the PFEM presented until now are the price to pay in order to gain all the benefits that this Lagrangian technique can give. The PFEM in fact makes simple some tasks that are extremely critical for a complex computational analysis and, in some cases, it allows the solution and the modeling of problems that other methods cannot even face.

The PFEM is designed for those problems where a huge change of topology occurs. The mesh nodes define the evolution of a domain during its history of deformation. In this way, the detection of the free surface in a fluid flow or the contours of a melting solid object do not introduce in the scheme any additional complication or implementation work because their contours are known at the beginning of each time step.

With the PFEM the bodies are totally free to deform and move in the space. There is not any limitation in this sense because the domain is continuously updated. Hence the final configuration can be very different then the initial one and occupy a volume of space arbitrary bigger than the initial one. This cannot been done with a classical Eulerian strategy, where a bounding box needs to be defined at the beginning of the analysis.

Furthermore with the PFEM the interface between a fluid stream and a solid is detected automatically via the same Alpha Shape method that generates the mesh. The interface is formed by the nodes of those hybrid elements (elements which nodes belong to both the solid and the fluid) that have been generated by the Delaunay triangulation and they fulfill the Alpha Shape criteria.

PFEM is also easy to couple with other numerical methods. This thesis is an example of the coupling of the PFEM with the FEM. In other works, it has been shown that the PFEM can also be easily coupled with discrete methods IEEEexample:PFEMbederosion. Furthermore, within the same PFEM, one may chose the preferred numerical method for solving the governing equations after the remeshing. In this work a mixed velocity-pressure stabilized formulation is used for the fluid, but one may implement any other Lagrangian methodology.

All this explains why in the past the PFEM has been employed for facing challenging simulations as tunneling IEEEexample:PFEMexcavation,IEEEexample:PFEMtunnelling, forming processes IEEEexample:PFEMforming,IEEEexample:nuestroPaperMan, melting of polymers IEEEexample:PFEMmelting,IEEEexample:PFEMmelting2, transport erosion and sedimentation in fluids IEEEexample:PFEMfsi, and fluid-multibody interaction IEEEexample:PFEMbedErosion2.

In conclusion, the PFEM is an extremely powerful technology which shortcomings are legitimized by the complexity of the problems to solve.

3.4.2 Mass conservation analysis

This section is focused on the analysis of the mass preservation feature of the PFEM-FIC stabilized formulation for free surface flow problems.

It is useful to distinguish two different sources of variation of volume (in all the examples presented in this work, the density is a constant, so the concepts of mass variation and volume variation are totally equivalent). The first source of volume variation is induced by the numerical solution for each time step increment. The volume may vary during the non-linear iterations due to the inaccuracy of the scheme. So the volume of the fluid domain changes from a volume computed at the beginning of the time step to a slightly different volume computed after the convergence of the solution scheme. This volume variation caused by the computation is called .

The second source of mass loss is related to the remeshing. After the achievement of the convergence the fluid domain is remeshed. During this process the volume may vary because some elements are erased or new elements are added to the mesh. The volume at the end of the time step is and the variation induced by the remeshing with respect to is named . In Figure 38 a graphical representation of the scheme is given.
Monitoring of volume within a time step interval (ⁿt,ⁿ⁺¹t).
Figure 38: Monitoring of volume within a time step interval .

The sum of both contributions gives the total volume variation for each time step such as

(227)

where

(228)

(229)

For the proper understanding of the mass preservation feature of a formulation is essential to distinguish these different mass losses sources. Nevertheless, it is improper to define these two sources as uncorrelated. In fact, the variation of mass generated by the remeshing may affect the mass losses of the non-linear scheme. In fact, the remeshing induces a perturbation of the equilibrium reached in the previous time step through the elimination and creation of elements. This may affect the convergence of the numerical scheme and, consequently, the satisfaction of the incompressibility constraint.

The numerical examples presented in this section are taken from IEEEexample:nuestroPaper. All the examples have been run using the real value of the bulk modulus in matrix (Eq.(220)). This matrix is the component of the tangent matrix for linear momentum equations that takes into account the pressure variation. In Section 3.4.3 it will be shown that substituting in matrix (Eq.(220)) the real bulk modulus with a pseudo bulk modulus defined as , is helpful for the conditioning and the overall accuracy of the analysis.

3.4.2.1 Numerical examples

In order to verify the mass preservation feature of the formulation, several problems for which guaranteeing the mass conservation is critical, have been solved. In this section, some of these involving impact and mixing of free surface fluids, are presented.

The problem has been solved first in 2D. Figure 39 shows the geometry of the tank, the material properties, the time step size and the initial mesh of 5064 3-noded triangles discretizing the interior fluid. The fluid oscillates due to the hydrostatic forces induced by its original position. In this and in the next problems, the effect of the surrounding air has not been taken into account.
Error creating thumbnail: File missing
Figure 39: 2D analysis of sloshing of water in a prismatic tank. Initial geometry, analysis data and mesh of 5064 3-noded triangles discretizing the water in the tank.

This problem, as for all the ones presented in this section, has been solved using the real bulk modulus in matrix (Eq.(220)), that is equivalent to compute the pseudo bulk modulus with the parameter .

Figures 40 and 41 show some snapshots of the water stream at different times. Pressure contours are plotted over the deformed configuration.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 40: 2D sloshing of water in a prismatic tank. Snapshots of water geometry at two different times (). Colours indicate pressure contours (I/II).
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 41: 2D sloshing of water in a prismatic tank. Snapshots of water geometry at two different times (). Colours indicate pressure contours (II/II).

The accumulated volume loss (in percentage versus the initial volume) for the method proposed is approximately 1.33% over 20 seconds of simulation time (curve of Figure 42). This value has been computed by summing all the volume variations due to the computation () for all the time steps of the analysis as

(230)

where is the initial volume.

In Figure 42 also the mass losses obtained using a standard first order fractional step method and the PFEM IEEEexample:nuestroPaper are plotted. Clearly the method proposed in this work leads to a reduced overall fluid volume loss.
Error creating thumbnail: File missing
Figure 42: 2D sloshing of water in a prismatic tank. Time evolution of the percentage of water volume loss due to the numerical algorithm. Comparison with the fractional step solution.

In some cases, it may be required to ensure the absolute absence of volume variations. With a small shifting of the free surface it is possible to guarantee the perfect mass conservation in the fluid domain. The method consists on moving the free surface nodes in the normal direction with an offset equal to , where the variation of volume is computed as the sum of the volume variation due to the remeshing at the previous time step and the volume variation due to computation of the current time step . is the length of the free surface. In Figure 43 a graphical representation of the technique in 2D is given.

Error creating thumbnail: File missing
Figure 43: Strategy for recovering the initial fluid volume by correcting the free surface at mesh generation level. (a) Compute total volume variation () before remeshing. (b) Compute free surface offset = . (c) Move free surface nodes in the normal direction to the boundary a distance equal to the offset computed in (b).

This procedure induces an arbitrary alteration of the equilibrium configuration, thus it has to be used cautiously. A fundamental requirement for using this technique is that the solution scheme produces very small volume variations and so the correction is minimal.

Figure 44 shows that the total fluid volume loss can be reduced to almost zero by applying this strategy at each time step.
Error creating thumbnail: File missing
Figure 44: 2D sloshing of water in a prismatic tank (). Mass losses applying the strategy described in Figure 43.

Figure 45 shows a comparison of the fluid volume losses between using the real bulk modulus () and an arbitrarily reduced pseudo bulk modulus (). For both cases, the time step increment is .

Error creating thumbnail: File missing
Figure 45: 2D sloshing of water in a prismatic tank. Time evolution of percentage of water volume loss obtained using the current method with (curve A) and (curve B) .

Figure 46 shows that a similar improvement in the volume preservation can be obtained using and reducing the time step to . This, however, increases the cost of the computations.

Error creating thumbnail: File missing
Figure 46: 2D sloshing of water in a prismatic tank. Time evolution of percentage of mass loss obtained with the current method. Curve A: and . Curve B: and .

These results show that accurate numerical results with reduced volume losses can be obtained by appropriately adjusting the parameter in the tangent bulk modulus matrix, while keeping the time step size to competitive values in terms of CPU cost. In Section 3.4.3 this issue will be analyzed in detail. Specifically, the effects of the bulk modulus in terms of volume preservation and overall accuracy will be analyzed and an efficient strategy to predict the best value for the pseudo bulk modulus () will be presented.

In Figure 47 the input data for the 3D case are given. The same sloshing problem of the 2D simulation is solved using a relative coarse initial mesh of 106771 4-noded tetrahedra and .
Error creating thumbnail: File missing
Figure 47: 3D analysis of sloshing of water in prismatic tank (). Initial geometry and analysis data.

Figure 48 shows the results for the 3D analysis at the same time instants of Figure 40.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 48: 3D analysis of sloshing of water in prismatic tank (). Snapshots of water geometry at two different times (s left and s right).

In Figures 49 and 50 the mass losses are analyzed. The values of Figure 49 have been computed as Eq.(230). It is remarkable that the percentage of total fluid volume loss due to the numerical scheme after 10 seconds of analysis is approximately 1%

The values plotted in the graph of Figure 50 are the percentage mass loss due to computation for each time step.

These values have been computed for each time step as

(231)
Error creating thumbnail: File missing
Figure 49: 3D analysis of sloshing of water in prismatic tank (). Time evolution of accumulated water volume loss (in percentage) due to the numerical algorithm.
Error creating thumbnail: File missing
Figure 50: 3D analysis of sloshing of water in prismatic tank (). Volume loss per time step over of analysis. Average volume loss per time step: 1.64%

This problem simulates the 2D motion, impact and subsequent mixing of two fluid streams originated by the collapse of two water columns located at the end sides of a prismatic tank. Figure 51 shows the initial geometry and the problem data.
Error creating thumbnail: File missing
Figure 51: Collapse and impact of two water columns in a prismatic tank. Analysis data, initial geometry and discretization (3988 3-noded triangles).

The two water columns have been discretized with 3988 3-noded triangles.

The effect of the surrounding air has not been taken into account in the analysis. The problem has been solved for .

Figure 52 shows four snapshots of the motion of the water columns after removal of the retaining walls. Alter a few instants the two water streams impact with each other and mix as shown in the figures.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 52: Collapse and impact of two water columns. Snapshots of the evolution of the flow at different times.
The evolution of the percentage of the initial fluid volume loss over the simulation time is shown in Figure 53.
Error creating thumbnail: File missing
Figure 53: Collapse and impact of two water columns. Accumulated fluid volume (in percentage) over eight seconds of analysis due to the numerical algorithm. Results for .

The values plotted in Figure 53 have been computed using Eq.(230). A maximum of 2.8% of the initial fluid volume is lost over eight seconds of analysis. This can be considered a low value for a problem of this complexity.

In Figure 54 the mass losses for each time step computed with Eq.(231) are plotted.
Error creating thumbnail: File missing
Figure 54: Collapse and impact of two water columns. Volume loss (in %) per time step. Results for .

The final example is the 3D analysis of the impact and mixing of a water drop as it falls in a cylindrical tank filled with the same fluid. Figure 55 shows the material and analysis data and the initial discretization (88892 4-noded tetrahedra).
Error creating thumbnail: File missing
Figure 55: Falling of a water sphere in a tank filled with water. Analysis data, initial geometry and discretization (88892 4-noded tetrahedra).

The problem has been solved using the real bulk modulus also in the tangent matrix ( taking ). Figure 56 shows four snapshots of the mixing process at different times.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 56: Falling of a water sphere in tank containing water. Evolution of the impact and mixing of the two liquids at different times. Results for .

The total water mass lost in the sphere and the tank due to the numerical algorithm was 2% after 3 seconds of analysis (Figure 57). The average volume loss per time step was .

Error creating thumbnail: File missing
Figure 57: Falling of a water sphere in a tank containing water (). Accumulated volume over three seconds of analysis due to the numerical algorithm.

3.4.3 Analysis of the conditioning of the solution scheme

This part is devoted to study the effect of the bulk modulus in the iterative matrix (matrix (Eq.(205)) on the partitioned solution scheme for quasi-incompressible fluids. For clarity purposes, the decomposition of matrix is recalled. is computed as

(232)

with

(233)

(234)

(235)

(236)

In this section, it is shown that matrix (Eq.(236)) can deteriorate the conditioning of the linear system (Eq.(204)) and the overall accuracy of the numerical scheme if the real bulk modulus is used. In order to avoid these drawbacks, a pseudo bulk modulus should be used. In this section, a practical rule to set up the value of a pseudo-bulk modulus a priori in matrix for improving the conditioning of the linear system is presented. The efficiency of the proposed strategy is tested in several problems analyzing the advantage of the modified bulk tangent matrix for the stability of the pressure field, the convergence rate and the computational speed of the analyses. The technique has been tested on the FIC/PFEM Lagrangian formulation presented, but it can be easily extended to other quasi-incompressible stabilized finite element formulations. The method proposed is based on using a pseudo-bulk modulus in the volumetric component (Eq.(236)) of the linear momentum tangent matrix (Eq.(232)), while the actual physical value of the bulk modulus is used for the numerical solution of the mass conservation equation, namely in matrices (Eq.(198)) and (Eq.(199)). The pseudo-bulk modulus is defined “a priori” as a proportion of the actual bulk modulus of the fluid (i.e. with ). The study here presented recalls the ideas presented in IEEEexample:nuestroPaperBulk.

This section is structured as follows. In the first part, the numerical inconveniences induced by the real bulk modulus are analyzed. Then the strategy for predicting the optimum value for pseudo bulk modulus is explained. Finally several numerical examples are given in order to validate the technique.

3.4.3.1 Drawbacks associated to the real bulk modulus

The linearized system (Eq.(204)) suffers from numerical instabilities due to the ill-conditioning of the iteration matrix (Eq.(232)) caused by the presence of the bulk modulus in matrix (Eq.(236)), so in the part of that takes into account the pressure variation. This problem was already pointed out in previous works where similar partitioned schemes were used IEEEexample:PavelPaperB,IEEEexample:PavelPaper,IEEEexample:PavelPaperA. The ill-conditioning of the iterative matrix of the linear momentum equations originates from the different orders of magnitude of its two main contributions: the mass matrix (Eq.(234)) and the bulk matrix (the contributions of the viscous matrix (Eq.(235)) and the geometric part (Eq.233)) are negligible, for low viscous flows). Typically the terms of the bulk matrix are orders of magnitude larger than those of the mass matrix.

A reliable measure of the quality of a matrix is the IEEEexample:BatheLibro. For a general matrix , the condition number is defined as

(237)

where denotes here the L2 norm of matrix .

The condition number gives an indication of the accuracy of the results from the matrix inversion and the linear equation solution. Values of close to 1 indicate a well-conditioned matrix.

The deterioration of the quality of matrix affects directly the convergence of the iterative linear solver. In this work, the iterative (BCG) solver has been used and its tolerance has been fixed to .

The time step increment affects highly the conditioning of the iterative matrix , as it has been shown in the first numerical example of Section 3.4.2.1. In fact, is inversely proportional to the time increment, while depends linearly on it (see Eqs.(234) and (236)). For this reason, is well-conditioned only for a tight range of time increments.

For the sake of clarity, a numerical example is used to visualize and quantify the inconveniences caused by the ill-conditioned matrix in the iterative solution scheme. The problem chosen is the 2D water sloshing in a prismatic tank presented in Section 3.4.2.1. However, with the purpose of reducing the computational cost of the analysis, the problem has been solved using a coarse mesh. So the initial geometry and the material data are the same as in Figure 39, but the mesh is the one shown in Figure 58a. In Table 58 all problem data are summarized.
Error creating thumbnail: File missing
(a) 2D water sloshing. Initial geometry and finite element mesh.
Figure 58: 2D water sloshing. Problem data.

To highlight the importance of the time step on the conditioning of the iterative matrix, the problem has been solved for two different time increments without reducing the modulus in matrix (). Using a time step of =, the iterative matrix has a condition number while, for = , the condition number is . If a larger time step is used, the linear system cannot even be solved. This deterioration is reflected by the number of iterations of the linear BCG solver: for the former case the average number of iterations is around , while for the latter it is around . Clearly, this also leads to a significant increase of computational time for solving a time step.

The ill-conditioning of the linear system also affects the rate of convergence of the iterative loop of the scheme given in Box 8.

In the graph of Figure 59 the convergence rates for the pressure and the velocity fields at =1.75 are displayed.
Error creating thumbnail: File missing
Figure 59: 2D water sloshing (=1). Convergence of the velocities and pressure at t=1.75 s.

For each iteration of the non-linear loop, the convergence rates for the nodal velocities and pressures have been computed as

(238)

(239)

where and are prescribed error norms for the nodal velocities and the nodal pressures, respectively. In the examples solved in this work, has been chosen.

The values of Figure 59 have been obtained considering a time step increment of =. The curves show that the convergence of the scheme is slow, especially for the pressure field. In fact, the pressure error (Eq.(239)) after 25 iterations is still larger than the pre-defined tolerance of . The lack of a good convergence for the pressure field has two main negative effects: the pressure solution is not accurate and mass conservation is not preserved.

The pressure contours at are shown in Figure 60.

Error creating thumbnail: File missing
Figure 60: 2D water sloshing (=1). Pressure contours at .

Concerning the mass conservation of the fluid, Figure 61 shows the accumulated mass variation in absolute value versus time. This values have been computed at each time step as

(240)

where is the volume loss due to computation and is the initial volume.

After 20 seconds of simulation, the percentage of mass loss is , which corresponds to a mean volume variation of per time step.

Error creating thumbnail: File missing
Figure 61: 2D water sloshing (=1). Accumulated percentage of mass variation in absolute value versus time.

In conclusion, the use of real bulk modulus () in matrix leads to ill-conditioned linear systems and, as consequence, it limits the range of suitable time step increments and deteriorates the convergence of the solution scheme.

These drawbacks were overcome in previous works IEEEexample:PavelPaper,IEEEexample:PavelPaperA by substituting the physical bulk modulus with a smaller pseudo bulk modulus . In this way, it is possible to extend the applicability of the partitioned scheme to a larger range of time increments. However, this strategy is based on heuristic criteria and cannot be used widely because each problem requires a specific value for . In other words, a pseudo bulk modulus that works well for certain analysis can fail for a different one. In this work, the optimum pseudo bulk modulus is computed and the strategy for predicting its value will be explained in the next section.

3.4.3.2 Optimum value for the pseudo bulk modulus

In this work, the pseudo bulk modulus is computed as choosing for the parameter a value able to guarantee the well-conditioning of the tangent matrix (Eq.(236)). The numerical examples presented in this section will show that this reduced value for the volumetric compressibility improves the overall accuracy of the scheme and does not affect the mass conservation. The reason is that the modification only affects the quality of the iterative matrix and the rate of convergence of the linear momentum equations, while the continuity equation, where the mass conservation constraint is imposed, is not modified. In other words, the physical value of the bulk modulus is used in matrices (Eq.(198)) and (Eq.(199)). This feature represents an innovation versus previous approaches where the pseudo bulk modulus is also used in the mass conservation equation IEEEexample:PavelPaperB,IEEEexample:PavelPaper.

The parameter is computed as

(241)

where the operator () denotes the mean of the absolute values of the non-zero matrix components.

For a uniform mesh of elements of characteristic size , the here called 'optimum value' of is estimated as follows

(242)

where is the value of the shape function at the element center and the following approximation has been used:

(243)

If water is considered ( and ) and linear triangles are used (), has the following dependency with the mesh size and the time step,

(244)

Typically, the parameter is calculated at the beginning of the time step at the first iteration of the non-linear loop. It can also be computed at every time step, at certain instants of the analysis or only once at the beginning of the analysis. This is so because the order of magnitude of does not vary during the analysis, unless the time step is changed, or a refinement of the mesh is performed. In these cases, needs to be calculated again because it has a square dependency on both parameters and . The numerical results presented in this work have been obtained computing via Eq.(241) at the beginning of the analyses only.

The parameter can be computed and assigned locally to each element or globally to the whole mesh. If it is calculated globally, all the elemental bulk contributions to the iteration matrix will have the same value of . Otherwise Eq.(241) is used considering separately each element of the mesh; each element will have a different value of . The former approach has a reduced computational cost but it works worse than the local approach for not uniform finite element discretizations. In particular, it is recommended to use the local approach when a refinement of the mesh is performed (in the next section, a problem with a refined zone is studied). In this work, unless otherwise mentioned, the global approach for computing is used.

3.4.3.3 Numerical examples

In order to show the efficiency of the method described in Section 3.4.3.2, two representative free surface problems involving the flow of water have been solved: the sloshing problem introduced in Section 3.4.3.1 and the collapse of a water column against a rigid obstacle presented in IEEEexample:FSImono1. The so-called dam break problem has been chosen here to demonstrate that the strategy does not affect the incompressibility constraint at all and the method is able to simulate problems of impact of fluids. For both sloshing and dam break problems, a comparison of the performances of the scaled bulk matrix and the scheme with is given.

The applicability and generality of the method is studied using very different average mesh sizes and time steps for both problems. Note that the dynamics of the sloshing problem is completely different from the dam break one.

The problem of Figure 58a is solved using the parameter computed as in Eq.(241). For , , while for =, . For both time step increments the resulting condition number of the iterative matrix is and the average number of iterations of the linear BCG solver is around 15. The improvement versus using is evident, if compared to the numbers presented in Section 3.4.3. Reducing the number of iterations of the linear solver, also reduces the computational time. For example, using for a duration of the sloshing simulation of , the total computational time for is while for it reduces to . Also the convergence of the non-linear loop improves with the optimum value of .

In Figures 62 and 63 the convergence of the velocity and pressure fields, respectively, obtained with and , is compared. The faster convergence of the solution using a smaller value of is noticeable.
Error creating thumbnail: File missing
Figure 62: 2D water sloshing. Convergence of the velocities at for and (optimum value).
Error creating thumbnail: File missing
Figure 63: 2D water sloshing. Convergence of pressure at for and (optimum value).

In Figure 64 the pressure solution at time t=1.75 obtained with =0.0053 is illustrated. It can be appreciated a remarkable enhancement versus the solution for =1 (see Figure 60). Note that the elements generated in the free surface region adjacent to the boundaries are due to the coarseness of the mesh and the remeshing criteria and not to the computation. A better result can be obtained using a smaller average mesh size or a refined mesh in the free surface region.

Error creating thumbnail: File missing
Figure 64: 2D water sloshing (). Pressure contours at .
As stated in Section 3.4.3, the convergence in the continuity equation affects the conservation of mass of the fluid. It has been just shown that the convergence of the pressure improves using a good prediction of . Consequently, also the mass conservation is better ensured using than using .
Error creating thumbnail: File missing
Figure 65: 2D water sloshing. Accumulated mass variation in absolute value along the duration of the analysis. Solution for and (optimum value).

In the graph of Figure 65 the percentage of accumulated mass variation (in absolute value) versus time obtained with and are compared. The values plotted in the graph have been computed according to Eq.(240). The better mass preservation of the solution with the smaller value of is remarkable.

In the graph of Figure 66 the accumulated mass variation obtained with is illustrated. The values plotted in the graph have been computed using Eq.(230). After 20 seconds of simulation, the scheme with the reduced value of has an accumulated mass variation of , which corresponds to a mean volume variation for each time step of . The solution with guarantees a better conservation of mass than with . This represents another evidence that the (quasi)-incompressibility constraint is not affected by reducing the bulk modulus for solving the momentum equations.

Error creating thumbnail: File missing
Figure 66: 2D water sloshing (=0.0053). Accumulated mass variation versus time.

In order to verify the applicability of the method, the sloshing problem has been solved using = and different mesh sizes. In particular, the following average mesh sizes have been used: =0.1, 0.15, 0.2, 0.3 and 0.4.

The problem was solved setting =1 and computing its reduced value using Eq.(241).

The curves in Figure 67 show that the reduced value of guarantees better results versus using . For all the meshes, the number of iterations required by the linear solver to reach a converged solution is smaller. Furthermore, the results show that the strategy is applicable to coarse and fine meshes.

Table 6 collects all the data and the results. The number of iterations of the linear solver has been considered as a quality indicator of the analyses. As shown in the previous sections this value is related to the condition number of the iterative matrix.

Error creating thumbnail: File missing
Figure 67: 2D water sloshing. Number of iterations of the linear solver for different numbers of velocity degrees of freedom. Results for =1 and the optimum value of .

Table. 6 2D water sloshing. Numerical values of the graph of Figure 67.
average degrees of freedom number of iterations
mesh size (velocities) =1 optimum






2D water sloshing. Numerical values of the graph of Figure 67.SL_meshes_iterationsT

The problem of Figure 58a has been solved using different time steps: 0.0001, 0.0005, 0.001, 0.005, 0.01, 0.02. The mesh has a mean size of 0.15. The numerical results obtained with the reduced value of and by setting =1 are compared. Once again, the number of iterations of the linear solver is the parameter chosen to indicate the quality of the analyses: the smaller this value is, the better conditioned the linear system is.

The graph of Figure 68 shows that the accuracy of the method does not depend on the time step increments, when the suitable reduced value for is used.
Error creating thumbnail: File missing
Figure 68: 2D water sloshing. Number of iterations of the linear solver for different time step increments. Results for =1 and the optimum value of .

Table 7 summarizes the problem data and collects the numerical values of the graph of Figure 68.

Table. 7 2D water sloshing. Numerical values of the graph of Figure 68.
number of iterations
=1 optimum







2D water sloshing. Numerical values of the graph of Figure 68.SL_coarse_stepsT

For each value of , the number of iterations is 16 and the condition number does not change. Furthermore, using a reduced value of allows us to solve the problem for each time increment, while if is fixed to 1, the results are acceptable only until =0.005. For larger time steps the results for =1 are not accurate or the analyses do not even converge.

As mentioned in the previous sections, the parameter has a square dependency on the mesh size (see Eq.(242)). For this reason, if the discretization is not uniform the global estimation of might not guarantee the well-conditioning of the linear system. In these cases, a local computation of is recommended. The sloshing problem is here solved again using the mesh shown in Figure 69 and . The spatial discretization has two different mean sizes: in the center =0.1 while =0.4 elsewhere.
Error creating thumbnail: File missing
Figure 69: 2D water sloshing. Finite element mesh with a refined zone.

The problem is solved both using =1 and by computing its optimum value globally and locally. The mean number of iterations required by the linear solver to converge for each of these three options is 319, 21 and 17 respectively. Hence, the local approach guarantees the best result for this type of non-uniform meshes.

In Figure 70 the solution at time obtained with the local approach is illustrated.

Error creating thumbnail: File missing
Figure 70: 2D sloshing with a refined zone. Solution at time obtained with the local approach.

In this section, the collapse of the water column induced by the instant removal of a vertical wall is studied.

Error creating thumbnail: File missing
(a) 2D dam break. Initial geometry.
Figure 71: 2D dam break. Problem data.

The initial geometry of the problem is illustrated in Figure 71a. In Table 71 all the problem data are collected. As for the previous example, the problem is first solved with a very coarse mesh. Then a comparison with the solution obtained for =1 is given. After that, the same problem is solved varying the mean mesh size and for different time step increments. The results obtained with the optimum value of are compared to the experimental results presented in IEEEexample:DBtest. The objective is to show that the reduction of the bulk modulus in the iteration matrix does not affect the numerical solution of this class of impact problems which can be solved also with a larger time step.

The problem is first solved with a coarse discretization (mean element size =0.0125). The solutions obtained with the optimum value of and with =1 are compared in terms of the condition number of the iteration matrix and the number of iterations of the linear solver. For =s, matrix has condition numbers =1028 and =60, for =1 and =0.0535, respectively, which correspond to 1251 and 14 iterations of the linear solver, respectively. For the same discretization and a time increment of =, Eq.(241) yields =0.000535. The condition number of and the iterations of the linear solver are the same as using a time step increment ten times smaller. Conversely, for =1, the condition number grows to =102090 and the iterative scheme does not converge.

The problem of Figure 71a has been solved for = using unstructured meshes with the following mean element sizes: = 0.004, 0.005, 00075, 0.01, 0.0125. The problem was solved both setting =1 and computing its optimum value using Eq.(241). As for the previous section, the number of iterations of the linear solver has been considered as the quality indicator of the analysis. In Figure 72 and Table 8 all the data and the results are collected.

Error creating thumbnail: File missing
Figure 72: 2D dam break. Number of iterations of the linear solver for different numbers of velocity degrees of freedom. Results for =1 and the optimum value of .

Table. 8 2D dam break. Numerical values of the graph of Figure 72.
average degrees of freedom number of iterations
mesh size (velocities) =1 optimum






2D dam break. Numerical values of the graph of Figure 72.DB_coarse_meshes

Figure 72 and Table 8 confirm that the efficiency of the method is not affected by the mesh size. In other words, the well-conditioning of the iterative matrix is guaranteed for coarse and fine meshes.

The same problem was solved for different time steps: =0.00005, 0.0001, 0.0005, 0.001, 0.0025. The mesh used has a mean size of =0.005. Table 9 collects the problem data and the resulting number of iterations for all the analyses solved with =1 and with the optimum value of .

Figure 73 and Table 9 confirm the conclusions of the previous section.
Error creating thumbnail: File missing
Figure 73: 2D dam break. Number of iterations of the linear solver for different time steps. Results for =1 and the optimum value of .

Table. 9 2D dam break. Numerical values for the graph of Figure 73.
numer of iterations
=1 optimum






2D dam break. Numerical values for the graph of Figure 73.DB_coarse_steps

The strategy is not affected by the time step and, contrary to the case for =1, it does not impose limitations to the range of time step increments for solving the problem. For example, for the present problem the time step increment is constrained by the geometry and the dynamics of the problem only. In other words, the maximum time step increment is the one that guarantees that the fluid particles do not cross through the boundaries. For the present problem and the chosen mesh, the maximum time step increment was =s.

All the examples presented in this section have shown the positive effect for the accuracy of the numerical scheme given by the use of the predicted value of in matrix (Eq.(232)). In all the following problems presented in this work, this strategy has been used. In particular, the technique is still available in FSI problems for improving the conditioning of fluid counterpart of the global tangent matrix.

3.5 Validation examples

This section is devoted to the validation of the unified stabilized formulation for fluids and solids at the incompressible limit. For each example several discretizations are analyzed in order to verify the convergence of the method. The formulation will be validated by comparing the numerical results to experimental tests and the numerical results of other formulations. In the first part quasi-incompressible Newtonian fluids are analyzed, then a problem involving a hypoelastic quasi-incompressible structure is studied.

3.5.1 Validation of the Unified formulation for Newtonian fluids

The purpose of this section is to validate the Unified formulation for the analysis of Newtonian free surface flows.

In all the examples presented in this section, in matrix (Eq.(236)) the pseudo bulk modulus is considered and its value has been computed according to the strategy presented in Section 3.4.3.2.

The effect of the boundary conditions is studied in detail. In particular, the strategy for modeling the slip conditions in the Lagrangian way described in Section 3.4.1.3 is tested and validated.

All the numerical examples have been solved for several meshes in order to verify the convergence of PFEM stabilized formulation.

In this section, the sloshing of a viscous fluid in a prismatic tank is analyzed. The initial configuration of the problem is illustrated in Figure 74a and the problem data are given in Table 74.

Error creating thumbnail: File missing
(a) Sloshing of a viscous fluid. Initial geometry of 2D simulation.
Figure 74: Sloshing of a viscous fluid. Problem data.
Slip conditions have been imposed on the tank walls. The problem has been solved in 2D with different discretizations in order to verify the convergence of the numerical scheme. In particular, the following average mesh sizes of 3-noded triangles have been used: 0.012m, 0.011m, 0.01m, 0.009m, 0.008m, 0.007m, 0.006m and 0.005m. In Figure 75 the finest (mesh size=0.005m) and the coarsest (mesh size=0.012m) meshes are illustrated.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
average mesh size= 0.012m
Figure 75: Sloshing of a viscous fluid. Coarsest (1967 triangles) and finest (11440 triangles) meshes used for the analysis.

All the analyses have been run for and the total duration of the study is .

In Figure 76 some representative snapshots of the 2D simulation are given. Each snapshot refers to the maximum height reached by the fluid during its sloshing. The pressure contours have been plotted over the deformed configuration. The pictures in Figure 76 assess the smoothness of the pressure field.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 76: Sloshing of a viscous fluid. Pressure contours over the deformed configuration at some time instants.
In Figure 77 the time evolution of the fluid level height at the left side of the tank is plotted for the two finest meshes used.
Error creating thumbnail: File missing
Figure 77: Sloshing of a viscous fluid. Time evolution of the free surface position at the left side of the tank. Solutions obtained with the finest meshes: average mesh size 0.006m and 0.005m.

The graph shows that a converged solution has been reached because the two diagrams almost coincide.

The convergence analysis has been performed for the time evolution of the free surface position at the left side of the tank. For the convergence study, the solution obtained with the finest mesh (average size 0.005m) has been considered as the reference solution. The error is computed as

(245)

where is the value computed for the finest mesh at the time step and is the total number of time steps.

For a given mesh and for each time step the fluid height at the left side of the tank is compared to the one obtained with the finest mesh at the same time instant. The graph of Figure 78 shows an overall quadratic slope, despite some oscillations in the curve.
Error creating thumbnail: File missing
Figure 78: Sloshing of a viscous fluid. Convergence analysis for the time evolution of the free surface level at the left side of the tank. Error computed with Eq.(245).

A reason for these oscillations is the local character of this convergence study. In order to reduce these drawbacks, a convergence study has been performed also for a global parameter, specifically the potential energy.

The potential energy has been computed as

(246)

where is the element mass and is the mean value for the nodal heights.

The graph of Figure 79 is the convergence curve for the potential energy error. Once again, the error is computed with respect to the results obtained with the finest mesh in the same form as in Eq.(245).
Error creating thumbnail: File missing
Figure 79: Sloshing of a viscous fluid. Convergence graph for the potential energy. Error computed according to Eq.(245).

Comparing this graph to the curve of Figure 78, the oscillations are reduced and the quadratic convergence for the used error measure is confirmed.

In this section, the collapse of a water column on a rigid horizontal plane is studied. The initial configuration of the problem is illustrated in Figure 80a and the problem data are given in Table 80.

Error creating thumbnail: File missing
(a) Collapse of a water column. Initial geometry for the 2D simulation.
Figure 80: Collapse of a water column. Problem data.
The phenomena involved in this problem are the collapse of a free surface water column and the consequent spread of the water stream over a horizontal plate. The numerical solution is compared to the experimental results given in IEEEexample:Martin where many experimental observations of the collapse of free liquid columns are collected. The characteristic variable chosen for the comparison is the residual height of the water column (see Figure 81).
Error creating thumbnail: File missing
Figure 81: Graphical definition of parameter .

The problem has been solved for both stick and slip conditions. The 2D simulation has been run for three different discretizations. In particular, the average mesh sizes , and have been used. In Figure 82 the finest (mesh size=) and the coarsest (mesh size=) mesh are illustrated.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
average mesh size= 0.006m
Figure 82: Collapse of a water column. Coarsest (258 triangles) and finest (13224 triangles) meshes used for the 2D analysis.

In 3D the problem has been solved only for a mesh of 338839 4-noded tetrahedra with average size 0.002m.

All the analyses have been run for .

In Figure 83 the velocity and the pressure fields obtained with the finest mesh and slip boundary conditions are plotted over the deformed configuration of the fluid.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Velocity field at Velocity field at
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Pressure field at Velocity field at
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Pressure field at
Figure 83: Collapse of a water column. Pressure and velocity contours over the deformed configuration at three time instants.
The results for the 3D simulation with slip conditions are shown in Figure 84.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Velocity and pressure fields at Velocity and pressure fields at
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Velocity and pressure fields at
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 84: Collapse of a water column. 3D simulation: velocity contours plotted over the 3D geometry and pressure contours drawn over the central section.

In Figure 85 the time evolution of the height at the left wall obtained with the finest mesh is compared to the values measured in the experimental tests IEEEexample:Martin.

Error creating thumbnail: File missing
Figure 85: Collapse of a water column: residual height vs time. Comparison between the best numerical solution and the experimental results.

The graph shows a very good agreement between the numerical and the experimental results.

For this problem, the solution obtained for slip and stick boundary conditions have been compared. In Figure 86 the velocity contours obtained for both types of boundary conditions are plotted over the deformed configuration for some time instants. The mesh used is the coarsest one (average size 0.006 m).

Error creating thumbnail: File missing
Error creating thumbnail: File missing
slip c., velocity field at slip c., velocity field at
Error creating thumbnail: File missing
Error creating thumbnail: File missing
stick c., velocity field at slip c., velocity field at
Error creating thumbnail: File missing
Error creating thumbnail: File missing
stick c., velocity field at
Figure 86: Collapse of a water column. 2D results for slip and stick conditions and a coarse mesh (average size 0.006m). Velocity contours over the deformed configuration at three time instants.
In Figure 87 the pressure field obtained for both boundary conditions is given.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
slip c., pressure field at slip c., pressure field at
Error creating thumbnail: File missing
Error creating thumbnail: File missing
stick c., pressure field at slip c., pressure field at
Error creating thumbnail: File missing
Error creating thumbnail: File missing
stick c., pressure field at
Figure 87: Collapse of a water column. 2D results for slip and stick conditions and a coarse mesh (average size 0.006m). Pressure contours over the deformed configuration at three time instants.

The stick condition affects a layer that has a size of the same order of magnitude than the discretization. Hence, the coarser the mesh is, the bigger is the zone affected by the stick condition. When a coarse mesh is used, as in this case, imposing the stick condition penalizes excessively the motion and it can provoke the instability illustrated in the right column of Figure 87. The stream obtained with the stick conditions is less uniform than the one obtained in the slip case and there are pressure concentrations. On the contrary, with a slip condition a fine solution is obtained also for the pressure field despite the coarseness of the mesh used in this example. This is because the motion of the wall particles helps in the remeshing step and in the global motion of the fluid stream.

In Figure 88 the results for the velocity and the pressure fields obtained with the mesh with average size and stick conditions are given. The pictures show that also for the stick case, a good solution for the pressure field can be obtained if the mesh is sufficiently fine.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
velocity field at velocity field at
Error creating thumbnail: File missing
Error creating thumbnail: File missing
pressure field at velocity field at
Error creating thumbnail: File missing
Error creating thumbnail: File missing
pressure field at
Figure 88: Collapse of a water column. Results obtained with stick condition and a mesh with average size of 0.00075m. Velocity and pressure contours over the deformed configuration at three time instants.
The graph plotted in Figure 89 shows that, also with the coarsest mesh tested in this example (average mesh size ), the solution obtained with the slip condition is very close to the experimental results, while the numerical solution obtained for the same mesh but imposing stick conditions on the walls is quite far from the expected one.
Error creating thumbnail: File missing
Figure 89: Collapse of a water column. Results for slip and stick conditions for the coarset tested mesh (average size ). Time evolution of the residual height of the water column.
However, as it is shown in the graphs of Figure 90, the solution obtained with the stick condition tends to the slip solution with mesh refinement. In particular, using a discretization with an average mesh size of the two solutions almost coincide.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
vs , mesh size= 0.00045m
Figure 90: Collapse of a water column. Results for slip and stick conditions for different average mesh sizes. Time evolution of the residual height of the water column.

The collapse of a water column over a rigid step is here studied in 2D and 3D. The initial configuration of the problem is illustrated in Figure 91a and the problem data are given in Table 91.

Error creating thumbnail: File missing
(a) Collapse of a water column over a rigid step. Initial geometry of 2D simulation.
Figure 91: Collapse of a water column over a rigid step. Problem data.

The phenomena involved in this problem are the collapse of the water column, the impact against the rigid step, the subsequent creation of the wave, the impact against the vertical wall and the final mixing of the fluid. The numerical solution is compared to the experimental results given in IEEEexample:DBtest, where many experimental observations of the collapse of free liquid columns are collected.

The problem has been solved for both stick and slip conditions. A convergence study for different mesh sizes has been performed for the 2D simulation. In particular, for the convergence test the following average mesh sizes of 3-noded triangles have been used: 0.011m, 0.0095m, 0.009m, 0.0085m, 0.008m, 0.0075m, 0.007m, 0.0065m, 0.006m, 0.0055m, 0.005m, 0.0045m, 0.004m, 0.003m, 0.0025m and 0.002m. In Figure 92 the finest (mesh size=0.0002m) and the coarsest (mesh size=0.0011m) mesh are illustrated.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
average mesh size= 0.011m
Figure 92: Collapse of a water column over a rigid step. Coarsest (1014 triangles) and finest (24668 triangles) meshes used for the 2D analysis.
All the analyses have been run for .
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 93: Collapse of a water column over a rigid step. Velocity contours over the deformed configuration at time: , , and . Results obtained with the finest mesh and slip conditions on the walls.

In Figure 93 the numerical results obtained for the 2D simulation with the finest mesh are compared to the experimental observations IEEEexample:DBtest for the same instants. The velocity contours are plotted over the deformed configurations.

The 3D problem has been solved considering a width of 0.146m and an average size of 0.008m for the tetrahedra of the FEM mesh. This gives an initial mesh of 170711 4-noded tetrahedra. In Figure 94 the 2D results (green) are superposed to the cut performed at the center of the 3D domain (grey) for the same time instants. An average size of 0.008m has been used for the edges of both the triangles and the tetrahedra for the 2D and 3D discretizations.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 94: Collapse of a water column over a rigid step. Superposition of the results obtained with an average mesh size of 0.008 m for the 2D and 3D simulations (green and grey colours, respectively) and with slip conditions on the walls.

The Figure 95 refers to the slip problem only and it shows that the solutions obtained with the three finest meshes are very similar. This is a further evidence of the convergence of the PFEM strategy here presented.

Error creating thumbnail: File missing
Figure 95: Collapse of a water column over a rigid step. Superposition of the 2D results obtained at with the meshes with average size 0.002m (orange), 0.0025m (pink) and 0.003m (green) (slip conditon on the walls).

In Figure 96 the solution obtained at with a stick condition and the finest mesh is given. The differences between the streams are very small.

Error creating thumbnail: File missing
Figure 96: Collapse of a water column over a rigid step. Superposition of the 2D results obtained at with the mesh with average size 0.002m for the slip (orange colour) and the stick (pink) conditions.

As already mentioned, the convergence analysis has been performed for both cases of slip and stick conditions imposed on the walls. Specifically, the study of convergence has been performed first for the maximum velocity and then the kinetic energy of the whole domain. Both values have been computed at , so just before the impact with the rigid step.

In the graph of Figure 97 the Y-axis is the maximum velocity obtained for the stick and slip conditions and the X-axis refers to the degrees of freedom () for the velocity. The solutions obtained with the stick condition suffer from oscillations for the coarsest meshes, but for finer meshes tend to those of the slip condition.
Error creating thumbnail: File missing
Figure 97: Collapse of a water column over a rigid step. Maximum velocity at for slip and stick conditions (2D simulation).

A convergence study for the kinetic energy of the whole domain at has been performed. The kinetic energy has been computed as

(247)

where is the element mass and is the mean value for the nodal velocities of the element.

The graph of Figure 98 is the evolution of the kinetic energy for the first 0.1 seconds of analysis, for both cases of slip and stick conditions and for three different meshes: the coarsest one (mesh size=0.011m), the finest one (mesh size=0.002m) and an intermediate discretization (mesh size=0.006m).
Error creating thumbnail: File missing
Figure 98: Collapse of a water column over a rigid step. Time evolution of the kinetic energy for slip and stick conditions and for three different FEM discretizations (2D simulation).

The curves for the slip case are almost superposed, contrary to the ones for the stick problem. However, for the finest mesh, the curves for the slip and stick conditions are very close.

In the graph of Figure 99 the Y-axis refers to the kinetic energy at for the stick and slip conditions and the X-axis refers to the degrees of freedom for the velocity. With slip conditions, the kinetic energy has almost the same value for all the discretizations (the difference between the values given by the coarsest and the finest discretization is lower than 0.5% However for stick conditions the kinetic energy grows with the refinement of the mesh and it seems to converge to the value of the slip case.
Error creating thumbnail: File missing
Figure 99: Collapse of a water column over a rigid step. Kinetic energy at for slip and stick conditions.
The convergence curve for the kinetic energy evolution on time for the first of analysis is given in Figure 100. For both slip and stick cases, the error has been computed using to Eq.(245). The slope is almost 2 for both curves.
Error creating thumbnail: File missing
Figure 100: Collapse of a water column over a rigid step. Convergence analysis for the kinetic energy evolution on time for the first of analysis for slip and stick conditions (2D case).

In Figure 101 the time evolution of the pressure until is studied. The graph on the left refers to a coarse mesh and the one in the right to the finest mesh. The plots show that the boundary conditions affect the pressure field. In particular the stick condition induces oscillations in the pressure solution which amplitude reduces by refining the mesh.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 101: Collapse of a water column over a rigid step. Time evolution of the pressure for a mesh with average size 0.006 m (top) and 0.002 m (bottom). Solution for stick and slip conditions (2D simulation).

This analysis has shown that for coarse meshes there are significant differences between the solution given by imposing stick or slip conditions. However, for fine meshes the problems solved with either slip or stick conditions tend to the same solution. In fact, for coarse meshes, the stick conditions affect highly the motion of the fluid and a refinement is required in order to obtain a solution close to the expected one. Instead, with slip conditions on the walls it is possible to obtain fine results, for both the velocity and pressure fields, also for coarse meshes.

3.5.2 Validation of the Unified formulation for quasi-incompressible hypoelastic solids

The Unified formulation is here validated for quasi-incompressible solids by solving a benchmark problem for non-linear solid mechanics, the Cook's membrane. In Section 2.3.4 the same structure was analyzed considering a compressible material. In this section, the nearly-incompressible problem is analyzed. In particular, first the problem is solved for with all the solid elements derived in this thesis, the V, the VP and the VPS elements comparing the numerical results for the displacements, the pressure and the stresses. Then the maximum displacement given by the three elements solving the same problem for the range of Poisson ratio from 0.3 to 0.499999 is compared.

In Figure 102a the initial geometry and the problem data are given.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
(a) Initial geometry
Figure 102: Nearly incompressible Cook's membrane. Initial geometry, material data and FEM mesh (127 elements).

The problem has been solved in 2D for different discretizations in order to verify the convergence of the scheme. In this case, unstructured meshes have been used. In Figure 102 the coarsest finite element discretization (average size 5) used in this example is depicted. The finest mesh tested in this problem has mean size 0.25 that gives a mesh of 52372 triangles.

As for the previous case, the problem is studied in statics and the self-weight of the membrane has not been taken into account. The reference solution for this problem is taken from IEEEexample:CMinc, where the problem was solved with an Incompatible Bubble element. In this publication the tip vertical displacement obtained by other formulations is also given. Specifically the results are provided for a FEM formulation with linear displacements and constant pressure and the Enhanced Assumed Strain formulation IEEEexample:CMinc2,IEEEexample:CMinc3. For all these formulations and the finest mesh tested in IEEEexample:CMinc, this value is around .

The Poisson ratio 0.4999 represents a material that is almost incompressible. It is well known that values of the Poisson ratio close to 0.5 generate numerical problems to non-stabilized FEM schemes, such as the locking of the solution. Apart from this, the proximity to the incompressible limit also produces ill-conditioned matrices and deteriorates the convergence of the solution. In order to overcome these drawbacks a stabilized formulation is required. In this work, the stabilization of the quasi incompressible mixed formulation for the solid element is achieved via the FIC strategy presented in Section 3.1 using the VPS-element. The quasi-incompressible membrane is here also solved using the non-stabilized V and VP elements derived in Chapter 2. Both these elements, although in different ways, suffer from instability near the incompressible limit of the material. Despite that, this example is presented with the purpose of showing how each formulation is affected by the mentioned drawbacks and to underline the superiority of a mixed formulation for dealing with nearly incompressible materials.

Error creating thumbnail: File missing
Figure 103: Nearly incompressible Cook's membrane. Top corner vertical displacement for the V, VP and the VPS elements.

In Figure 103 the top corner vertical displacements obtained with the V, the VP and the VPS elements for different FEM meshes are plotted.

In Table 10 the numerical values are given.

0.9

Table. 10 Nearly incompressible Cook's membrane. Top corner vertical displacement for different formulations and discretizations.
mesh number of V-element VP-element VPS-element
size elements
5 127 4.411 7.031 7.268
4 194 4.365 7.178 7.338
3 361 4.648 7.401 7.508
2 802 5.643 7.551 7.603
1.5 1441 4.937 7.632 6.655
1 3288 5.690 7.695 7.714
0.75 5806 6.199 7.707 7.729
0.5 13015 6.731 7.745 7.755
0.25 52372 7.272 7.765 7.771

Nearly incompressible Cook's membrane. Top corner vertical displacement for different formulations and discretizations. CMincTable

For the V, VP and VPS the percentage errors for the tip displacement with respect to the reference solution are 5.68%, 0.707 % and 0.791%, respectively.

For the analysis of stress results, the case of the mesh with average size 1.5 is studied. In Figure 104 the X-component of the Cauchy stress tensor obtained with V, VP, and VPS elements is plotted over the deformed configurations. From the pictures it is clear that both the V and VP solutions deteriorate for values of close to the incompressible limit of the material, while the stress field given by the VPS-element is good and smooth.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
V-element VPS-element
Nearly incompressible Cook's membrane. Results of the XX component of the Cauchy stress tensor for V, VP and VPS elements.
Figure 104: Nearly incompressible Cook's membrane. Results of the XX component of the Cauchy stress tensor for V, VP and VPS elements.
In Figure 105 the pressure solution obtained with the VP-element and the VPS-element are given.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
VP-element
Figure 105: Nearly incompressible Cook's membrane. Pressure solution for the VP and VPS elements.

The non-stabilized VP-element yields a pressure field that is completely untrustworthy exhibiting the classical checkerboard modes. Instead, the solution of the VPS element is smooth and accurate.

From the kinematic point of view, the displacements obtained by the mixed formulation (VP and VPS elements) are close to the expected solution, while the solution given by the V-element is totally locked. For the mesh displayed in Figures 104-105 (average mesh size 1.5) the errors for the top corner displacement with respect the reference solution are 36.5%, 1.9% and 1.61% for the V, the VP and the VPS elements, respectively.

The same problem is solved for different Poisson ratios, from 0.3 to 0.499999, using the same FEM discretization (average mesh size equal to 1). For the stabilized formulation the Poisson ratio that appears in the tangent matrix of the linear momentum equations has been limited to 0.4999. Instead, in the equation of the pressure the actual value of the Poisson ratio is used. This technique is similar to the one used for the stabilization of the fluid in Section 3.4.3 where the bulk modulus is reduced only in the tangent matrix of the linear momentum equations but not in the continuity equation. The graphs of Figure 106 are the plots of the top corner vertical displacement obtained for the different Poisson ratios using the V, the VP and the VPS elements (the graph at the right is a zoom of the last part of the graph at the left).
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 106: Nearly incompressible Cook's membrane. Maximum vertical displacement for different Poisson ratios. Results for V, VP and VPS elements. The graph at the right is a zoom of the last part of the graph at the left.

In Table 11 the numerical values are given.

0.9

Table. 11 Nearly incompressible Cook's membrane. Numerical values for the maximum vertical displacement for different Poisson ratios obtained with V, VP and VPS elements.
Poisson ratio V-element VP-element VPS-element
0.3 9.176 9.212 9.245
0.4 8.557 8.601 8.628
0.45 8.143 8.203 8.230
0.49 7.658 7.833 7.866
0.495 7.503 7.781 7.817
0.499 6.988 7.735 7.766
0.4999 5.690 7.695 7.714
0.49999 4.872 7.532 7.714
0.499999 4.706 6.696 7.713

Nearly incompressible Cook's membrane. Numerical values for the maximum vertical displacement for different Poisson ratios obtained with V, VP and VPS elements. CMpoiTable

The graphs show that for Poisson ratios bigger than 0.45, the velocity formulation exhibits locking. Instead, the mixed formulations despite the increasing of ill-conditioning of the linear system, yieldsa good solution until a Poisson ratio of 0.49999. Beyond this value, the non-linear solver does not even converge for the non-stabilized VP-element, while the solution given by the VPS-element is still a good one.

3.6 Summary and conclusions

In this chapter the Unified Stabilized formulation for quasi-incompressible materials has been derived. This numerical procedure is based on the mixed Velocity-Pressure formulation derived in Chapter 2 for compressible materials.

In order to deal with material incompressibility using a linear interpolation for both the pressure and the velocity fields, the scheme required to be stabilized. The necessary stabilization has been introduced using an enhanced version of the Finite Calculus (FIC) method. In Section 3.1 the complete derivation of the stabilized form of the mass balance equation has been presented. The FIC stabilization has been derived for quasi-incompressible fluids and has been extended also to hypoelastic quasi-incompressible solids. The solid finite element based on the mixed Velocity-Pressure stabilized formulation has been called VPS element.

The solution schemes for solving quasi-incompressible Newtonian fluids and hypoelastic solids with the stabilized mixed Velocity-Pressure formulation have been given and explained in detail.

Section 3.4 has been fully devoted to the analysis of free surface fluids with the Unified stabilized formulation. First the Particle Finite Element Method (PFEM) has been explained and its advantages and disadvantages have been highlighted. A simple technique for modeling the slip conditi

ons moving the wall particles have been also given. Then mass preservation properties of the PFEM-FIC stabilized formulation have been tested with several numerical examples. It has been shown that the method yields excellent results for a variety of 2D and 3D free surface flow problems involving surface waves, water splashing, violent impact of flows with containment walls and mixing of fluids. Next the conditioning of the scheme has been studied and the numerical inconveniences associated to the high value of the bulk modulus have been highlighted. An efficient and easy to implement technique for improving the conditioning and the global convergence of the problem using a scaled value for the pseudo bulk modulus has been given. The strategy has been validated for two benchmark problems for free surface flows.

In the last section of this chapter several numerical examples have been presented for validating the unified stabilized formulation for solving fluids and solids close to the incompressible limit. The proposed method has been validated versus both experimental tests and numerical results from other formulations and it has been shown that the method is convergent to the expected solution. Particular attention has been given to the analysis of the boundary conditions. In particular, it has been shown that for inviscid fluids and coarse meshes it is preferable to use slip conditions in order to avoid the pressure concentrations induced by the stick conditions.

4 Unified formulation for FSI problems

Chapter 4. Unified formulation for FSI problems

4.1 Introduction

In the previous chapters the Unified formulation was used for solving the mechanics of fluids and solids separately. In this chapter it will be shown that the method is also suitable for solving fluid-structure interaction (FSI) problems. The algorithm for coupling the three solid formulations (the V, the VP and the VPS elements) with the mixed Velocity-Pressure FIC-stabilized formulation for Newtonian fluids is described. It will be shown that the proposed method gives the possibility to select the formulation for the solid depending on the problem to solve. The numerical results will be validated by comparing the solution of benchmark FSI problems with the results of the literature.

In this work, the FSI problems are solved in a monolithic way. This means that fluids and solids are solved within the same linear system and the iterations between the fluid and the solid solver are not required, as for staggered schemes. The unified formulation for both fluid and solid mechanics solution uses the same framework (Lagrangian description), the same spatial discretization and the same temporal integration schemes (linear shape functions and implicit integration of time), the same unknown variables (velocities and pressures), and the same solution method (two-step partitioned scheme). In conclusion, solids and fluids just represent different regions of the same global domain and they differ on the characteristic material parameters only. All this simplifies the implementation of the FSI solver. In fact, it is not required any variable transformation and the implementation effort for coupling the mechanics of fluids and solids is reduced to a proper assembly of the global linear system.

Concerning the mesh, there must be a correspondence at the interface between the solid and the fluid nodes, in the spirit of conforming mesh methods. It will be shown that the PFEM guarantees the conformity.

In conclusion, for explaining the extension of the Unified formulation to FSI problems only the assembly of the global linear system and the way to detect the fluid-solid interface have to be described. This will be done in the next section. Then the solution schemes for solving FSI problems coupling the FIC-stabilized mixed VP formulation for quasi-incompressible Newtonian fluids with the V, the VP and the VPS elements for hypoelastic solids are given. Finally, several numerical examples of benchmark FSI problems are solved and all the schemes are validated and compared.

4.2 FSI algorithm

The assembly of the global linear system is performed by making a loop over all the nodes of the mesh. Each node provides the contributions of the elements that share the node, and each element is computed according to the specific constitutive law and solution scheme. So, when an interface node is analyzed, it is necessary to sum the contributions of both materials in the global linear system (see Figure 107 for a graphic representation of the global assembly).
Error creating thumbnail: File missing
Figure 107: Graphic representation of domain contributions to the global linear system.

Those nodes that belong to a single domain assemble to the global matrix with the elemental contributions of one material only.

Concerning the degrees of freedom, each node of the mesh is characterized by a single set of kinematic variables. So they move according to a unique velocity. This means that the degrees of freedom for the solid and fluid velocities coincide also at the interface nodes. On the other hand, in order to guarantee the correct boundary conditions for the stresses, each interface node has a degree of freedom for the pressure of the fluid and, for the VP and the VPS elements, for the pressure of the solid. This means that, for each interface node the fluid elements assemble only the contributions for the degree of freedom of the fluid pressure, while the solid elements do that for the solid pressure. This requires solving twice the continuity equation: once for the fluid domain and once for the VP-element or the VPS-element for the solid domain.

In order to ensure the coupling, the fluid and the solid meshes must have in common the nodes along the interface. In other words, there must be a node to node conformity. In this work, the solid is solved using the FEM while the fluid is computed by the PFEM. In terms of meshing, this means that the solid domain maintains the same grid during all the analysis, whereas the fluid is remeshed whenever its discretization becomes excessively distorted. The conformity of the fluid and solid meshes on the interface is guaranteed by exploiting the capability of the PFEM for detecting the boundaries. In practice, the fluid detects the solid interface nodes (nodes that are located on the external boundaries of the solid fixed mesh) in the same way it recognizes its rigid contours. As it has been explained in Section 3.4.1.1, this is done by an efficient combination between the Alpha Shape method and the Delaunay triangulation. According to this strategy, if the separation of the fluid contour from the solid domain is small enough so that that the Alpha Shape criteria are fulfilled, a fluid element connecting the fluid domain to the solid domain is generated. Otherwise the two domains keep apart from each other. In Figure 108 a graphic representation of the method for detecting the interface is given IEEEexample:Cremonesi.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
(a)
Figure 108: Detection of an interface with the PFEM IEEEexample:Cremonesi.

In the situation described in the pictures of Figure 108a none of the contact elements generated by the Delaunay triangulation fulfill the Alpha Shape criteria. Hence for the forthcoming time step interval there is not interaction between the solid and the fluid domains. Instead in Figure 108b some hybrid elements that share solid and fluid nodes have been generated. In this case, the coupling is active and the interface nodes assemble as described before.

For the good use of this procedure, it is essential to select a similar mesh size for the fluid and the solid domain in the interface zone. Otherwise, some typical drawbacks may arise. For example, if the solid mesh is much finer than the fluid one some distorted elements may form in the interface zone. The opposite case is even worse. If the solid mesh at the boundaries is much coarser than the fluid one, some fluid particles may cross the interface and compromise the whole computation. Using the same mesh size (and a reasonable time step increment) neither of these problems arise. Note that this behavior is characteristic of PFEM. In fact all this can occur also in the contact zone of a fluid domain with a rigid boundary in a fluid dynamic analysis, as explained in the Section 3.4.1. This confirms again that, from the remeshing point of view, there is not difference between the contour nodes of a deformable body and the nodes forming a rigid wall.

4.3 Coupling with the Velocity formulation for the solid

The coupling of the PFEM-FIC stabilized VP formulation with the V-element for the solid consists essentially on assembling adequately the linear momentum equations. For each degree of freedom for the interface nodal velocity, the solid and the fluid contributions sum. On the other hand, the stabilized continuity equation is assembled and solved only for the fluid elements because the pressure is not an unknown variable of the Velocity formulation.

In order to avoid ambiguities in the notation, the variables and the matrices referred to the fluid domain will be marked by subscript f ', and those ones belonging to the solid by s '. When a shared variable is considered (for example, the unknowns referred to the nodes of the interface) the subscript will be s,f '. Instead, the subscript that marks all the nodes of the domain is s+f '.

In Box 12, the solution scheme for a generic time interval is given. In the scheme only the contributions of the interface nodes are analyzed. For the rest of nodes the problem is locally uncoupled, hence the schemes described in Box 3 for hypoelastic solids and in Box 8 for Newtonian fluids still hold.

All the matrices and the vectors that appear in Box 12 were already defined in Boxes 4 and 9.

4.4 Coupling with the mixed Velocity-Pressure formulation for the solid

The coupling of the fluid stabilized VP formulation with the VP-element is performed similarly to the V-element. In particular, the solution scheme of the linear momentum equations does not change, because also in this case the degrees of freedom of the velocity are not duplicated. Hence the fluid velocity coincides with the solid velocity.

The main differences concern the pressure solution. The solid and the fluid pressures are two different degrees of freedom. As a consequence, the continuity equations are solved separately for the fluid and the solid. For the solid, both the stabilized and the not stabilized equations, Eq.(86) and Eq.(112) respectively, can be solved. Depending on the compressibility of the solid bodies one may choose one of the two schemes.

For the fluid counterpart, matrix (Eq.(220)) is computed using the pseudo bulk modulus and its value is predicted using the strategy described in Section 3.4.3.2.

In Box 13, the solution scheme for a FSI problem solved using the VP-element or the VPS-element for the solid is given for a time interval . Once again, all the matrices and vectors that appear Box 13 refer only to the contributions of the interface nodes and they have already been defined in Box 6 for the VP-element, Box 11 for the VPS-element and Box 9 for the stabilized VP formulation for quasi-incompressible Newtonian fluids.

For the nodes that do not belong to the interface, the schemes described in Box 5 and 8 still hold.

4.5 Numerical examples

This example is the two-dimensional abstraction of the moving of a circular cylinder between two parallel walls. The cylinder moves perpendicularly to its axis due to the gravity force increasing the falling velocity until an asymptotic value.

The distance from the rigid walls and the axis of the cylinder is . The radius of the circle is . The geometry of the problem and the material data are given in Figure 109a and in Table 109.
Error creating thumbnail: File missing
(a) Falling of a cylinder in a viscous fluid. Initial geometry.
Figure 109: Falling of a cylinder in a viscous fluid. Problem data.

The same problem has been studied in many other works IEEEexample:GilISPM,IEEEexample:GilISPM, IEEEexample:ImmersedWang. The numerical solution can be also compared to the analytical study of the motion of a rigid cylinder with constant velocity between two parallel plane walls. For this problem, the analytical solution is obtained studying the creeping motion equations using stream functions IEEEexample:HappelBook. Considering stick conditions on the walls, for a unit of length of the cylinder the drag force in stationary conditions is

(248)

where is the velocity of the rigid cylinder.

This relation for the drag force holds also for the stationary conditions of the problem studied in this section. From the equilibrium of the drag force with the Archimedes' force yields

(249)

where is the difference between the density of the two materials involved. Combining Eq.(248) and Eq.(249), the terminal velocity of a rigid circular cylinder with radius is

(250)

For this problem this relation gives

(251)

This value for the terminal velocity is near to the asymptotic ones obtained by the mentioned works IEEEexample:GilISPM,IEEEexample:GilISPM, IEEEexample:ImmersedWang.

Numerical studies show that this expression holds only for a tight range of values of the viscosity and, specifically, is not suitable for a fluid with small viscosity.

In this work, the solid has been modeled not as a rigid body but with a hypoelastic model, and using a high value for the Young modulus. The VP formulation has been used for the solid. This means that the example belongs properly to the class of the FSI problems. The problem has been solved for both stick and slip conditions on the vertical walls of the container. Due to the dominant role of the viscosity, the solution of the stick problem does not tend to the solution of the slip problem.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
t = 0.375 s t = 1.000 s
Error creating thumbnail: File missing
Figure 110: Falling of a cylinder in a viscous fluid. Snapshots of the cylinder motion at different instants of the 2D simulation (slip conditions on the boundaries). Velocity contours are depicted over the solid and fluid domains.
In Figures 110 and 111 some snapshots of the simulation with the velocity contours are given. Figure 110 refers to the slip case, Figure 111 to the stick one.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
t = 0.375 s t = 1.000 s
Error creating thumbnail: File missing
Figure 111: Falling of a cylinder in a viscous fluid. Snapshots of the cylinder motion at different instants of the 2D simulation (stick conditions on the boundaries). Velocity contours are depicted over the solid and fluid domains.

The numerical results plotted in Figures 110 and 111 have been obtained with a mesh with average size . The pictures show the importance of the boundary conditions for this case. It is evident that the velocity fields obtained for the slip and the stick cases, as well the terminal velocity of the cylinder, are significantly different.

For the same mesh, the resulting pressure field for the slip and stick cases is illustrated in Figure 112. The pictures show that the perturbation caused by the motion of the cylinder is almost imperceptible and there are not significant differences between the slip and stick cases.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
slip conditions, t=1 s
Figure 112: Falling of a cylinder in a viscous fluid. Pressure field obtained for the slip and stick cases.
In the graph of Figure 113 the time evolution of the vertical velocity of the cylinder obtained with the finest meshes (average size=) are given for both the slip and stick cases.
Error creating thumbnail: File missing
Figure 113: Falling of a cylinder in a viscous fluid. Time evolution of the vertical velocity of the cylinder. Results for the slip and the stick cases.

The terminal velocities of the cylinder obtained with slip and stick conditions are and , respectively.

For this example the transmission conditions between the solid and the fluid domain have been monitored. The curves of Figure 114 represent the time evolution of the Neumann condition in the X-direction (horizontal) at the points located at the boundary of the cylinder and depicted in Figure 109a. Specifically, the value plotted in the curves is the mean value of the X-component of vector () computed for the fluid and the solid elements at the points of Figure 109a.
Error creating thumbnail: File missing
Figure 114: Falling of a cylinder in a viscous fluid. Time evolution of the X-component of computed at the point of Figure 109a.

The graph shows that the transmission condition is guaranteed during all the analysis.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Mesh size= 0.0015m
Figure 115: Falling of a cylinder in a viscous fluid. Coarsest (average mesh size=0.0015) and finest (average mesh size=0.0004) meshes used for the example depicted on the solid domain.

A convergence study was also performed. In Figure 115 the coarsest (average mesh size=0.0015m) and the finest (average mesh size=0.0004m) meshes used for the example are depicted on the solid body. The whole domain has been discretized using 2983 and 42873 3-noded triangular elements, respectively.

The graph of Figure 116 shows the time evolution of the vertical velocity obtained with five of the meshes tested for the slip case.
Error creating thumbnail: File missing
Figure 116: Falling of a cylinder in a viscous fluid. Time evolution of the vertical velocity of the cylinder for five different meshes and slip conditions on the walls.
In Figure 117 the results obtained with the same meshes assuming stick conditions on the walls are given.
Error creating thumbnail: File missing
Figure 117: Falling of a cylinder in a viscous fluid. Time evolution of the vertical velocity of the cylinder for five different meshes and stick conditions on the walls.
The convergence graphs for both the stick and slip conditions are plotted in Figure 118. Both cases show the same convergence rate.
Error creating thumbnail: File missing
Figure 118: Falling of a cylinder in a viscous fluid. Convergence curves for the slip and stick cases.

The problem has been solved also for a quasi-incompressible solid using the VPS-element. For this case, a Poisson ratio of 0.4999 and the same Young modulus used for the problem solved with the VP-element, have been considered. In Figure 119 the velocity and the pressure fields for the solid and the fluid computed at for stick conditions on the walls and using a mean mesh size of are given.

In the graph of Figure 120 the time evolution of the vertical velocity obtained with the VPS-element for is compared with the solution obtained with the VP-element for and the same average mesh size and boundary conditions.

As expected, the solutions are almost the same and the vicinity to the incompressible limit does not compromise the quality of the results.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Velocity field Solid Cauchy stress
Error creating thumbnail: File missing
Figure 119: Falling of a quasi incompressible cylinder in a viscous fluid. Velocity field, fluid pressure and solid Cauchy stress (YY component) at for stick conditions on the walls and a mean mesh size of .
Error creating thumbnail: File missing
Figure 120: Falling of a cylinder in a viscous fluid. Solutions obtained using for the solid the VP-element () and the VPS-element ().

The problem was presented by Aristoff in IEEEexample:aristoff. In the mentioned work the experimental results of the water entry of spheres of different materials are studied.

In this section, the case of a nylon sphere is analyzed. The numerical results given by the Unified formulation (with the VP-element for the solid) are compared to the results of the laboratory test. The sphere impacts the water in the tank with a vertical velocity of 2.17 and its diameter is 2.54 cm. The density of nylon is 1140 and the Young modulus and Poisson ratio are 3 and 0.2, respectively. The water was modeled considering a density of 1000 , a dynamic viscosity 0.00089 and a bulk modulus 2.15 . In order to simulate this problem correctly, a very fine mesh is necessary. For this reason the whole domain was discretized with 1059924 tetrahedra. The time step used for the analysis is .

In Figure 121 the numerical results are compared to the experimental ones IEEEexample:aristoff. The comparison shows the good agreement of the results obtained with the PFEM Unified formulation with the laboratory values and its capability to simulate the complex phenomena such as the creation of the void (the air has not been taken into account) within the fluid domain.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 121: Water entry of a nylon sphere: comparison between experimental and numerical results.

This example is inspired from a similar problem presented in IEEEexample:Cremonesi. A volume of a viscous fluid drops from a rigid container over a thin and highly deformable membrane. The impact of the fluid mass causes an initial huge stretching of the structure and its subsequent oscillations. A hypoelastic law is used for the material of the structure. The problem was solved in 2D for two different values of the fluid viscosity, namely 50 and 100 , and with both the V and VP elements. The purpose was to compare the formulations and to show that both solid elements can be used for the modeling of standard elastic solids in FSI problems. The initial geometry of the problem is given in Figure 122a and the material data are given in Table 122a.
Error creating thumbnail: File missing
(a) Filling of an elastic container with a viscous fluid. Initial geometry.
Figure 122: Filling of an elastic container with a viscous fluid. Problem data.

In the graph of Figure 123 the results for the less viscous case (=50 Pas) obtained using the V and the VP elements for the solid are given. The comparison is performed for the vertical displacement of the lowest point of the elastic structure. The curves are almost coincident and only after of simulation some slight differences appear.

For the same problem, some representative snapshots are collected in Figure 124. The pressure contours are depicted over the solid domain, while over the fluid one the mesh is plotted. The numerical results refer to the simulation using the VP-element for the solid.

Error creating thumbnail: File missing
Figure 123: Filling of an elastic container with a viscous fluid (). Vertical displacement of the bottom of the elastic container obtained using the V and the VP elements for the solid domain.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 124: Filling of an elastic container with a viscous fluid (). Snapshots at different instants of the 2D simulation. Pressure contours depicted over the solid domain.

In Figure 125 the snapshots of the most viscous case (=100 Pas) are given for the same time instants of Figure 124. The numerical results refer to the solution obtained using the VP-element for the solid domain.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 125: Filling of an elastic container with a viscous fluid (). Snapshots at different instants of the 2D simulation. Pressure contours depicted over the solid domain.

Also for the most viscous case, the results obtained using the velocity and the mixed velocity-pressure formulations for the solid are compared analyzing the time evolution of the vertical displacement at the bottom of the elastic structure. The two curves of Figure 126 represent the solutions obtained with the velocity and the mixed velocity-pressure formulations.

Once again, the differences between the results of the two formulations are very small and this is a further evidence of the validity and flexibility of the proposed Unified formulation. This example evidenced the possibility to choose for the solid a velocity or a mixed formulation also for solving FSI problems.

Error creating thumbnail: File missing
Figure 126: Filling of an elastic container with a viscous fluid (). Vertical displacement of the bottom of the elastic container obtained using the V and the VP elements for the solid domain.

The problem illustrated in Figure 127 was introduced by Walhorn IEEEexample:walhorn.

Error creating thumbnail: File missing
Figure 127: Collapse of a water column on a deformable membrane. Initial geometry and problem data.

The water column collapses by instantaneously removing the vertical wall. This originates the flow of water within the tank, the formation of a jet after the water stream hits the rigid ground, and the subsequent sloshing of the fluid as it impacts a highly deformable elastic membrane. The membrane bends and starts oscillating under the effect of its inertial forces and the impact with the water stream.

The problem has been solved in 2D as in 3D, considering stick conditions for the rigid walls. In Figure 128 some representative snapshots of the 2D simulation are given.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 128: Collapse of a water column on a deformable membrane. Snapshots of the 2D simulation at different instants. The VP element is used for the solid.
The results obtained with the present formulation have been compared to the ones computed in [IEEEexample:FSImono1,IEEEexample:Idelshon,IEEEexample:Cremonesi]. In the graph of Figure 129 the time evolution of the horizontal deflection of the left top corner is illustrated.
Error creating thumbnail: File missing
Figure 129: Collapse of a water column on a deformable membrane. Horizontal deflection of the left top corner on time. Numerical results obtained with V and VP elements for the solid. Comparison with numerical results obtained in [IEEEexample:walhorn,IEEEexample:Idelshon,IEEEexample:Cremonesi].

The diagram shows that, for the first part of the analysis, the proposed formulation agrees well with the results reported in the literature. After around of simulation, the numerical results of each formulation starts to diverge, although for all the formulations the membrane oscillates two times around its vertical position before the time instant . The first part of the simulation is easier to analyze than the second one because the phenomena to model are less aleatory and the fluid splashes do not affect the results, as it occurs in the second part of the simulation. Furthermore, the initial deformation of the elastic structure affects highly the rest of the simulation. In fact, a smaller bending of the membrane induces an impact of the water stream at a higher height of the containing wall and with a bigger tangential component of the impact velocity. Consequently, the fluid stream impacts against the right side of the elastic membrane later and with reduced inertial forces. These considerations have an even greater effect for the 3D simulation.

In the graph of Figure 130 the results of the 2D and 3D simulations obtained using the V-element are compared.
Error creating thumbnail: File missing
Figure 130: Collapse of a water column on a deformable membrane. Horizontal deflection of the left top corner on time. Comparison between 2D and 3D analyses (V-element for the solid part).

The graph shows that there is a good agreement between the 2D and 3D analyses in the first part of the graph. Then the graphs are quite different. In particular, in the 3D problem the fluid stream impacts later against the right side of the membrane producing a reduced displacement of the elastic structure. This occurs because, in the 3D simulation, the fluid stream after the impact against the vertical rigid wall (this occurs after around of analysis), can move in the direction transversal to the main direction of its motion and this reduces the impact force of the water stream. This behavior is intensified if stick boundary conditions are considered and the mesh is not sufficiently fine.

In Figure 131 the numerical results of the 3D simulation are given.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 131: Collapse of a water column on a deformable membrane. Snapshots of the 3D simulation of different instants. The V-element is used for the solid.

4.6 Summary and conclusions

In this chapter the Unified Stabilized formulation has been used for the solution of FSI problems. Specifically, the numerical schemes for coupling the stabilized velocity-pressure formulation for quasi-incompressible Newtonian fluids with the V, the VP and the VPS elements for the solid have been presented. It has been shown that, depending on the specific need, one may choose any of these solid elements.

It has been shown that the implementation effort for extending the Unified formulation to the coupled FSI problems is small. Essentially, it consists of assembling adequately the global linear system and exploiting the capability of the PFEM for the detection of the interface.

Several validating examples have been given. The numerical solution obtained with the proposed scheme has been compared to analytical solutions and numerical results of other formulations for free surface FSI problems. Good agreement has been found in all cases. It has also been shown that the method is convergent.

5 Coupled thermal-mechanical formulation

Chapter 5. Coupled thermal-mechanical formulation

5.1 Introduction

This chapter is devoted to the coupling of the unified formulation with the heat transfer problem. The objective is to solve general problems involving fluids, solids or both of them which behavior is temperature-dependent and where the phase change of the materials is allowed.

There is a growing interest in industry in having a technology capable to solve coupled thermo-mechanical problems. In fact, industrial problems generally lack of an analytical solution and often the unique way to predict their solution is based on experience. An alternative is represented by experimental tests. However, these cannot be carried out for all the required problems; for example due to the complexity of the geometry, and, in many cases, the high cost associated to these. For these reasons, computer-based methods represent often the most reliable alternative for the analysis of coupled thermal-mechanical problems.

Thermo-mechanical problems involving fluid and solids however, represent a challenge for the numerical analysis due to its multi-field nature and high non-linearity. From the computational point of view, there are many complications associated to these problems.

First of all, the coupling with the heat problem increases the non-linearity of the mechanical problem. The solution of the mechanical problem depends on the heat transfer via temperature dependence of the material properties, at the same time, the heat problem depends on the solution of the mechanical problem due to the change of the configuration.

The complexity of the problem increases if phase change is also considered. This phenomenon not only increases the non-linearity of the problem due to the huge changes of topology that it may induce, but it also requires the implementation of a specific technology for its simulation. Furthermore accounting for phase change may induce problems related to the quality of the mesh. Consider, for example, the melting of a solid body. If a Lagrangian mesh is used, as it generally occurs with solids, the body that melts undergoes large deformations that may compromise the quality of the mesh and the results of the overall simulation.

Another critical point is the imposition of the boundary conditions for the temperature. For this it is required to track the contours of the domains involved in the simulation for all the duration of the analysis. This task may not be trivial if the problem involves a complex geometry, severe changes of topology or phase change.

With the PFEM many of these complexities are overcome thanks to its Lagrangian nature. In fact, as it has already explained in the previous chapters, the boundaries of the bodies are automatically detected by the position of the boundary nodes. Furthermore, if PFEM is also used for modeling the solid dynamics, the drawbacks associated to the large deformations of the solid caused by heat are overcome through the remeshing which guarantees a fine discretization for all the duration of the analysis. Finally, in the Lagrangian formulation the convective term disappears from the heat transfer equation. For all these reasons, other authors faced this problem in the past using the PFEM IEEEexample:PFEMmelting, IEEEexample:PFEMmelting2, IEEEexample:PavelThesis.

In this work, the PFEM is only used for the fluid domain. So in order to overcome the decay of quality of the solid mesh in a melting process, a further technology is required. The algorithm that has been used in this work consists on transforming the solid elements into fluid elements when they reach a threshold value. This means that the part of the solid body that is melting is considered as a fluid. Consequently it benefits from the PFEM technology and it can be remeshed if required. The details of the algorithm will be given in the following sections.

Due to the complexity and the vastness of the issue, in this treatise various assumptions and simplifications have been accepted. This has been in some cases at the expense of the accuracy of the analyses. Nevertheless, the main objective of this part of this work is to show that the PFEM unified formulation for FSI problems can be easily coupled to the heat problem in order to solve simulations where also the temperature affects. The principal assumptions of this work are described in the following.

In all the simulations, the air has not been taken into account, so heat transmission only occurs via the contact between solid and liquid bodies, or through the contours where the boundary conditions for the temperature have been imposed. Furthermore, in the phase change modeling, there is not a transition where a multi-fluid analysis is employed. The material that changes its state takes instantaneously the properties of the other material involved in the analysis, without considering a layer formed by a transition material. Note that these simplifications have been assumed just for convenience and they are not required by the formulation itself. With an additional implementation work both assumptions can be avoided. A proof of this are the successful attempts made in previous works where the PFEM was successfully used for modeling multi-fluids IEEEexample:PFEMmultiF.

In order to guarantee a strong thermal-mechanical coupling, the heat problem is solved within the same iteration loop used for solving the mechanical problem. This strategy increases the non-linearity of the mechanical problem because of its dependence on the temperature. A simpler way to solve the global problem is to consider a constant temperature during the non-linear iterations of the mechanical problem and to solve the heat problem only once the mechanical problem has converged. This approach leads to a lower computational cost and it reduces the non-linearity of the mechanical problem. However, the resulting coupling is weaker. Both strong and weak coupling strategies belong to the staggered solution schemes because the mechanical and the heat problem are solved in different linear systems. Conversely, in monolithic approaches both problems are solved simultaneously. On the one hand, this scheme enables the stability and the convergence of the whole coupled problem, on the other, it leads to longer and worse conditioned algebraic systems. However, in this work this scheme has not been chosen essentially for another reason. Monolithic schemes require modifications also in the formulation for the mechanical solution. This occurs because in the assembly of the linear system also the degrees of freedom of the temperature have to be taken into account. Instead, in staggered schemes, the FSI formulation remains the same because the mechanical and the heat transfer problems are two different blocks and they may be implemented separately. This, apart from reducing the implementation effort, gives us the possibility to test the unified formulation in its original version.

This chapter is split into the following sections. First, the governing equations of the heat are introduced and discretized using FEM. Also the scheme for the solution of a general time step is given. Then the thermal coupling algorithm is described and validated with three numerical examples. After that, the strategy for modeling phase change is explained. The chapter ends with an example where the thermal-mechanical coupling technique is applied.

5.2 Heat problem

The heat transfer problem in a Lagrangian setting is governed by the following differential equation

(252)

where is the temperature, is the thermal capacity, is the heat conductivity and is the heat source.

Notice that the convective term does not appear in Eq.(252).

The temperature and the heat flux at the boundaries are prescribed with the following boundary conditions

(253)
(254)

where and are the prescribed temperature and the prescribed normal heat flux at the boundaries and , respectively, and is the direction normal to the boundary.

The problem is completed by the initial condition

(255)

The space for the test functions for the temperature is defined as

(256)

Multiplying Eq.(252) by the test functions and integrating over the updated configuration domain, the following global integral form is obtained

(257)

Integrating Eq.(257) by parts, using Eq.(254) and neglecting the space changes of the conductivity, the following weak variational form of the heat problem is obtained

(258)

5.2.1 FEM discretization and solution for a time step

For the heat problem, the same procedure for the spatial discretization of the mechanical one is used. Hence the analysis domain is discretized into finite elements with nodes in the standard manner, leading to a mesh with elements and nodes. For 2D problems 3-noded linear triangles are used, while 3D domains are discretized using 4-noded tetrahedra. The temperature is interpolated over the mesh in terms of its nodal values, in the same manner as for the velocities and the pressure, using the global linear shape functions spanning over the elements sharing node (). In matrix form,

(259)

where and

(260)

In Eq.(259) vector contains the nodal temperatures for the whole mesh. Substituting Eq.(259) into Eq.(258) and choosing a Galerkin formulation with leads to the following algebraic equation

(261)

The matrices and vectors in Eq.(261) are assembled from the element contributions given in Box 1.

Eq.(261) is solved using a standard forward Euler scheme. For the iteration within a time interval the following linear system is solved

(262)

with

(263)

If the following conditions are verified, the next time step is considered.

(264)

where is a prescribed error norm.

5.3 Thermal coupling

In a thermo-coupled mechanical problem, the deformation induced by the heat has to be taken into account. The total deformation rate is the sum of the elastic, the plastic and thermal parts. Hence Eq.(122) for a thermal-elastoplastic problem is:

(265)

where the thermal part of the deformation rate tensor is defined as

(266)

where is the thermal expansion coefficient and .

In this work, the heat problem is solved after the mechanical problem in the same iteration loop (Figure 132).
Inner staggered thermal coupling scheme for a general time step (ⁿt,ⁿ⁺¹t).
Figure 132: Inner staggered thermal coupling scheme for a general time step .

This strategy belongs to the class of staggered schemes because the heat and the mechanical problems are solved in two different linear systems. This method is called ' scheme' in order to distinguish it from another staggered scheme that will be presented later. This strategy increases the non-linearity of the problem because the properties of the material may depend on the temperature. The convergence is tested at the end of each iteration for the velocities, the pressure and the temperature.

The thermal coupling can be also performed through a different staggered scheme, that will be called ' scheme'. In this case the heat problem is solved after the convergence of the mechanical problem (Figure 133).

External staggered thermal coupling scheme for a general time step (ⁿt,ⁿ⁺¹t).
Figure 133: External staggered thermal coupling scheme for a general time step .

The external staggered scheme reduces the non-linearity of the mechanical problem because during the non-linear iterations the temperature is considered as a constant. It also gives the possibility to choose different time step increments for the mechanical and the heat problems. On the other hand, the resulting coupling is weaker than for the internal scheme. For guaranteeing a stronger coupling the mechanical and the thermal problems should be solved iteratively within the time step. However, this strategy increases the computational cost of the analyses. In practice, their duration increases by a factor equal to the number of the global iterations of the thermal-mechanical problem.

5.3.1 Numerical examples

In this section, three numerical examples are presented. First, the heat transfer problem is validated comparing the numerical solution of a heated plate to the analytical one. Next, the sloshing of a fluid in a heated tank is presented. Finally, the heating of three solid objects in a fluid contained in a rectangular tank is studied.

The objective of this first example is to verify the implementation of the thermal problem. The problem consists of a square plate at initial temperature T=100 cooled by three of its edges at constant temperature T=0 keeping the other edge insulated. Figure 134 shows a graphical representation of the problem.

Error creating thumbnail: File missing
Figure 134: Heating of a plate. Initial geometry, thermal boundary and initial conditions.

The thermal conductivity of the plate is . The length of the plate edges is .

The problem can be computed analytically by solving the partial differential equations of the heat problem with the proper initial and boundary conditions. The analytical solution for the temperature field is:

(267)
In the graph of Figure 135 the temperature evolution with time at the central point of the top edge (with coordinates ) obtained with the proposed approach is compared to the analytical solution. Apart from the initial instants, the curves are almost identical.
Error creating thumbnail: File missing
Figure 135: Heating of a plate. Temperature evolution on time at the point . Analytical and numerical solutions.

This fluid dynamics problem was presented in IEEEexample:CoupledThermal. A fluid at initial temperature oscillates due to the hydrostatic forces induced by its initial position in a rectangular tank heated to a uniform and constant temperature of . The geometry and the problem data of the 2D simulation are shown in Figure 136.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 136: 2D sloshing of a fluid in a heated tank. Initial geometry, problem data, thermal boundary and initial conditions.

The fluid domain has been initially discretized with 2828 3-noded triangles. The coupled thermal-fluid dynamics simulation has been run for using a time step increment of =.

The main purpose of this example is to show the capability of the proposed Lagrangian technique for dealing with thermal coupled problems involving severe changes of topology. For this reason, some simplifications have been accepted. For example, the properties of the fluids do not depend on temperature. Furthermore, with the purpose of reducing the computational time and visualizing better the temperature contours, a very high (and not realistic) thermal conductivity has been used. These assumptions reduce clearly the truthfulness of the problem but, as it has already pointed out, this is not the principal objective of this example.

In Figure 137 the evolution of the temperature with time at the points and of Figure 136 is plotted. The coordinates of these sampling points are (), () and (), respectively.
Error creating thumbnail: File missing
Figure 137: 2D sloshing of a fluid in a heated tank. Time evolution of the temperature at the points and of Figure 136.

The fluid, because of its high thermal conductivity, changes its temperature quickly due to the contact with the hotter tank walls. The heat flux along the free surface has been considered to be null.

Figure 138 shows some snapshots of the numerical simulation. The temperature contours have been superposed on the fluid domain at different time instants.
Error creating thumbnail: File with dimensions greater than 12.5 MP
Figure 138: 2D sloshing of a fluid in a heated tank. Snapshots of fluid geometry at six different time instants. Colours indicate temperature contours.

Figures 138 and 137 show that the fluid does not heat uniformly because of the convection effect that is automatically captured by the Lagrangian technique here presented.

This 2D example is taken from IEEEexample:nuestroPaperMan. Three solid objects with the same shape fall from the same height into a tank containing a fluid at rest. For the solid objects a hypoelastic model has been used. The geometry and the problem data, as well the initial thermal conditions, are shown in Figure 139.
Error creating thumbnail: File with dimensions greater than 12.5 MP
Draft Nogué 457856805-ballDataGeo.png
Error creating thumbnail: File missing
Figure 139: Falling of three objects in a heated tank filled with fluid. Initial geometry, thermal conditions and material properties.

Once again, the material properties are assumed to be temperature independent. The fluid in the tank has initial temperature T=340, while the solid bodies from left to right, have initial temperatures T=180, T=200, T=220, respectively. The solid and the fluid domains have been discretized with a mesh composed of 9394 3-noded triangular elements. The simulation has been run for a total duration of using a time step increment =. The heat flux in the normal direction is assumed to be null for the boundaries in contact with the air or the walls.

In the graph of Figure 140 the evolution of the temperature at the central point of the three objects is plotted.
Error creating thumbnail: File missing
Figure 140: Falling of three objects in a heated tank filled with fluid. Evolution of the temperature at the center of the three objects.

Figure 141 collects six representative snapshots of the numerical simulation with the temperature results plotted over the fluid and the solid domains.

t=0.67st=1.82s t=2.66s
t=0.67s t=2.66s
t=4.50s t=7.00s
t=4.50s t=7.00s
t=8.00s Falling of three objects in a heated tank filled with fluid. Snapshots with temperature contours at different time steps.
t=8.00s
Figure 141: Falling of three objects in a heated tank filled with fluid. Snapshots with temperature contours at different time steps.

5.4 Phase change

For dealing with the phase change transformation, a term that takes into account the latent heat released or absorbed during the melting process need to be considered IEEEexample:NargesThesis. Hence, the heat equation for phase change problems reads

(268)
where is the latent heat of the phase transformation and is the liquid fraction function which is equal to one in the liquid phase and null in the solid one, as shown in Figure 142.
Error creating thumbnail: File missing
Figure 142: General liquid fraction function for phase change modeling.

The phase change is modeled in the following way.

The solid elements transform into fluid elements if at least one of the two following conditions is verified:

  1. : if the mean temperature of the solid elements that share the same node is greater than the melting temperature of the material. In this case, all these solid elements shift to the fluid domain, as illustrated in Figure 143 for a 2D problem.
  2. Error creating thumbnail: File missing
    Figure 143: Graphic representation of the change of phase algorithm.
  3. : if a solid element has accumulated a plastic deformation greater than a pre-defined limit value.

The transformation from fluid to solid or viceversa entails just a few of changes and small computational work. The main ones concern the change of the material properties and the remeshing procedure.

With respect to the former point, in this work the properties of the melted or solidified material changes instantaneously taking the values of the other material analyzed in the problem. This represents a strong simplification and it may affect the accuracy of the numerical results. For a better modeling of the problem, a multi-fluid analysis should be performed considering different properties for the melted or solidified materials. Despite this more accurate strategy has not been implemented or tested, from the author point of view, it can be easily linked to the unified formulation. This opinion is sustained by the fact that the proposed strategy is open to analyze different materials at the same time and the PFEM has already been used successfully in the past for solving multi-fluid problems IEEEexample:PFEMmultiF.

Regarding remeshing, as already explained in the previous chapters, the mesh is regenerated just on the fluid domain. So when a solid element passes to the fluid domain, its nodes, except the ones that still belong to the contour of the remaining part of solid domain, become involved in the remeshing procedure, as the rest of the elements in the fluid domain. The remaining part of the solid domain is not remeshed but it just loses the part of its domain involved in the melting.

The above two changes are essentially the only ones required by the thermally coupled unified formulation for dealing with phase change problems.

5.4.1 Numerical example: melting of an ice block

This 2D example is presented to show an application of the phase change strategy implemented. The problem was presented in IEEEexample:nuestroPaperMan. An ice block at initial temperature T=270 is dropped into a tank containing water at rest at temperature T=340. The initial geometry, the problem data and the initial thermal conditions are shown in Figure 144.

Error creating thumbnail: File with dimensions greater than 12.5 MP
Draft Nogué 457856805-iceDataGeo.png
Error creating thumbnail: File missing
Figure 144: Melting of an ice block. Geometry, material data and initial thermal conditions.

Ice is treated as a hypoelastic solid until some of its elements reach the fusion temperature (T=273.15). These elements pass to the fluid domain taking all its physical properties. For this analysis the following assumptions have been made: the mechanical and thermal properties of the water and the ice do not change with the temperature and the heat normal flux along the boundaries in contact with the air or the walls have been considered to be null.

In Figure 145 snapshots of some representative instants of the analysis are shown. The temperature contours are plotted over the water and the ice.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 145: Melting of an ice block. Snapshots with temperature contours at different time steps.
In Figure 146 the detail of the melting of the ice piece is illustrated. The finite element mesh is drawn over the solid and the fluid domains.
t=2.81st=4.70s t=9.77s
t=12.74s t=14.69s
t=15.47s Melting of an ice block. Zoom on the piece of ice at different time steps.
Figure 146: Melting of an ice block. Zoom on the piece of ice at different time steps.

5.5 Summary and conclusions

The purpose of this chapter has been to show that the Unified PFEM formulation can be easily coupled with the heat transfer problem in order to solve coupled thermal mechanical problems.

The coupling has been ensured via a staggered scheme for which the mechanical and the thermal problems are solved within the same iteration loop. For the temperature field the same linear shape functions of the velocity and pressure fields have been used.

Several numerical examples have been presented with the objective of showing the applicability of the Unified coupled thermal-mechanical formulation for solid and fluid dynamics problems involving the temperature.

The algorithm for the phase change modeling has been explained and an explicative numerical example has been also given.

The numerical examples presented have shown the possibility of the Unified PFEM formulation for dealing with thermal-mechanical problems.

6 Industrial application: PFEM Analysis Model of NPP Severe Accident

Chapter 6. Industrial application: PFEM Analysis Model of NPP Severe Accident

6.1 Introduction

In this chapter the analysis of an industrial problem solved using the Unified formulation is presented. The project required the modeling of the damages in a pressure vessel structure caused by the dropping of a volume of corium at high temperature.

In a nuclear power plant, the reactor pressure vessel (Figure 147) is the structure that contains the reactor core, the nuclear reactor coolant and the core shroud.

File:Draft Nogué 457856805-pressureVesselReal.png Pressure vessels for nuclear power plants. From IEEEexample:powerPlant.
Figure 147: Pressure vessels for nuclear power plants. From IEEEexample:powerPlant.

The core of a nuclear reactor is the place where the nuclear reactions occur. The coolant is used to remove the heat form the nuclear reactor and the core shroud is a strainless structure used to direct the cooling water flow.

Corium is a heterogeneous material, also called Fuel Containing Material (FCM) or Lava-like Fuel Containing Material (LFCM). Its formation is the result of the combustion and melting of the reactor's components and the subsequent chemical and radioactive reactions. Consequently, the corium's composition depends on the type of the reactor and specifically on the materials of its components. Generally it is composed by nuclear fuel, fission products, control rods, structural materials from the affected parts of the reactor, products of their chemical reaction with air, water and steam, and, in case the reactor vessel is breached, molten concrete from the floor of the reactor room.

The temperature of corium depends on its internal heat generation dynamics, specifically on the chemical and radioactive reactions that may occur during the accident. Potentially, the corium can reach temperatures over 2800 C.

Corium accumulates at the bottom of the reactor vessel and, in case of adequate cooling, solidifies. Otherwise, the corium may melt through the reactor vessel and flow out (Figure 148).
A blob of corium in the Chernobyl Nuclear Reactor. From IEEEexample:nuclearReactor.
Figure 148: A blob of corium in the Chernobyl Nuclear Reactor. From IEEEexample:nuclearReactor.

The seriousness of this kind of accident explains the high interest of this study.

The purpose of the analysis was to model with the PFEM the interaction between the corium and the pressure vessel structure. The scope was to evaluate the capability of the method for simulating such a complex problem. For the project, the PFEM unified formulation with thermal coupling presented in this work has been used.

The study is a multi-physics and highly non-linear problem and it involves many critical topics for a computational analysis, such as free-surface flow, fluid-structure interaction, plasticity with thermal softening and phase change.

Specifically, the analysis consisted of two models:

  • Basic model. A sphere of corium falls from a certain height over a prismatic plate;
  • Detailed model. A corium volume is placed all around a control rod housing located in the center of a steel shell.

The first model represents a general case and its geometry is not necessary related to a specific component of the pressure vessel. In Figure 149 a graphic and simplified representation of the study is given.

Basic model. Graphic representation of the phenomena required to model.
Figure 149: Basic model. Graphic representation of the phenomena required to model.

Instead, in the second model the geometry is more complex and it reproduces a part of the bottom head of a reactor pressure vessel. In Figure 150 a graphic representation of the area of study is shown.

Draft Nogué 457856805-pressureVessel.png
Error creating thumbnail: File missing
Figure 150: Graphic representation of the accumulation of corium at the bottom of the pressure vessel (images provided by NSSMC).

For both models, the phenomena to model are:

  1. Dropping and slumping of a high temperature volume of corium on the pressure vessel;
  2. Melting, deforming and distruction of the structure due to the heating caused by the corium;
  3. Slumping of corium through the vessel.

6.1.1 Assumptions allowed by the specification

Due to the complexity of the analysis the following assumptions were considered to be admissible in the specifications provided by the contractor for the analysis of the two models:

  • The data analysis given in the specifications (geometric data, material properties, boundary conditions...) were open to changes after consultation with the contractor;
  • Simplified constitutive models could be considered whenever necessary;
  • The dimensions of the models could be scaled down in consideration to the speed of the analyses;
  • The cooling water and its effects (boiling, flow etc.) are not considered in the study;
  • The corium is assumed to be a highly viscous fluid at constant temperature;
  • The corium could not mix or react with the surrounding media;
  • The melting and solidification of steel are treated with a simple model.

6.2 Numerical method

For simulating both models described in the previous section, the unified formulation with thermal-mechanical coupling has been used. The multi-physics of the project represented a great opportunity to verify the efficiency of the formulation in its completeness. In fact, the analysis involves most of the technical aspects presented in this work. The unified formulation allows the simulation of the interaction of free-surface fluids with solids with plasticity, as shown in Chapters 2, 3 and 4. Finally Chapter 5 showed that with just a small implementation work, also coupled thermal-mechanical and phase change problems can also be simulated.

From the practical point of view, in industrial projects the computational cost of the analysis has to be taken in serious consideration. Due to the complexity of the (3D) geometry, the duration of the physical phenomena and the required accuracy of the results, some industrial problems may lead to prohibitive computational times. Hence, in some cases it is necessary to accept a compromise between the accuracy and computational cost of the analysis. In order to calibrate at best the model, some preliminary studies of the simulations should be done. In this work, for example, this phase has been extremely useful for defining the adequate time step increment and the best FEM discretizations, deciding where a finer mesh was necessary and where not, and also for setting up some geometric and material values. For both models considered a preliminary study has been done and in the following some details of this phase are given.

Concerning the numerical method used for these simulations, the mechanical problem was solved with the unified formulation using the VP-element for the solid with an hypoelastic-plastic model, a von Mises yield criterion and thermal softening. The phase change was modeled following the criteria presented in Section 5.4. In both models there is only melting of the solids and not solidification of the fluids. For simplicity, it was assumed that the solid elements that become fluid take instantaneously the same material properties of the fluid. As already explained, for a more realistic simulation, a multi-fluid analysis should be performed taking into account different properties for the melted material.

6.3 Basic Model

6.3.1 Problem data

The basic model consists of a corium sphere at that falls from a height over a prismatic plate at .

The geometry of the problem is illustrated in Figure 151.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 151: Basic model: XZ and XY views of the initial geometry.

The boundary conditions of the plate are illustrated in Figure 152.

Error creating thumbnail: File missing
Figure 152: Basic model: constraints of the plate.

The properties of the corium are the ones given in the specification for and they are collected in Table 12. It has been assumed that the corium keeps the same temperature during the analysis.


Table. 12 Material data for the corium at T=2000K.
7865 22 500

Material data for the corium at T=2000K. MatDataCorium

Concerning the steel plate, the material data at are collected in Table 13. These material properties are comparable to those of the stainless steels and they can represent a Stainless Steel 309.


Table. 13 Material data for the plate at T=293 K.
7850 196 0.33 16 500 0.282 1693

Material data for the plate at T=293 K. MatDataPlate

The thermal influence on the material properties has been taken into account by multiplying the Young modulus , the heat conductivity , the specific heat and the yield stress by the coefficients , , and , respectively. Hence:

(269)
(270)
(271)
(272)

In Table 14 the values of the thermal coefficients are given for certain temperatures. For other temperatures, the coefficients are obtained via a linear interpolation.


Table. 14 Thermal coefficients for the plate.
Temperature Thermal coefficients
[C] [K]
20 293 1 1 1 1
93 366 0.992 1.082 1.036 1
204 477 0.969 1.208 1.091 1
315 588 0.914 1.333 1.146 1
427 700 0.910 1.460 1.202 1
538 811 0.864 1.586 1.257 0.781
649 922 0.828 1.711 1.312 0.767
760 1033 0.900 1.837 1.367 0.636
871 1144 0.778 1.962 1.422 0.267
982 1255 0.733 2.087 1.477 0.146
1093 1366 0.733 2.087 1.477 0.010
1377 1650 0.733 2.087 1.477 0

Thermal coefficients for the plate. TherDataPlate

6.3.2 Preliminary study

In order to calibrate the model, the problem was studied in 2D first. The original purpose of this preliminary study was just to understand better the difficulties of the problem and to determine a range of suitable values for both the mesh size and the time step increment. However, this preliminary study allowed to identify that some parameters of the proposed model needed to be modified in order to simulate the desired phenomenon of heating and melting of the steel plate.

In fact, in the initial configuration proposed in the specifications, the sphere of corium was located at the height . However, the preliminary studies showed that it would be not possible to simulate the desired problem if the corium would fall from that height. In fact, the inertial forces accumulated by the corium during its fall generate a huge splashing of the fluid over the plate. In order to allow the laying of the corium on the plate and the consequent heating of the structure, the initial height of the corium was reduced in the PFEM analyses. After some tests the height was chosen for the simulation.

The preparatory tests also showed that with the value of corium viscosity proposed in the specifications (), the corium flows progressively towards the edges of the plate and finally a part of this falls down by the plate contour. In order to avoid this undesired phenomenon, a non-linear viscosity was used for the basic model. In particular, corium was modeled using the constitutive law proposed and successfully applied for the simulation of melted metals in IEEEexample:ZincVisc2, IEEEexample:ZincVisc3. According to the mentioned publications, the viscosity of the corium was computed as:

(273)

where is the deviatoric strain invariant and is the yield stress.

6.3.3 Numerical results

The simulation was run using a finite element mesh composed by 327451 four-noded tetrahedral elements with 61452 nodes; 27936 nodes for the fluid and 33516 for the solid. The fluid was discretized using an uniform mesh with an average size of . On the other hand, the solid plate was not meshed uniformly. In its central part, for guaranteeing a good contact between the solid and the fluid, an unstructured mesh with the same average size of the fluid was used, while for the surrounding zone the average mesh size is .

For the first phase, when the splashing occurs the time step increment was used. In Figure 153 some snapshots referred to this phase are given.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 153: Basic model. Snapshots of the initial splashing of the corium on the plate.

In order to reduce the computational time of the analysis, for the other phases of the simulation the time step increment was increased to .

After the splashing and spreading of the corium on the steel plate, the plate starts heating due to the contact with the hotter corium and progressively melts in its central part. In Figure 154 a few snapshots of this phase are given.
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 154: Basic model. Snapshots of the melting of the plate and plot of temperature contours (blue and red colors correspond to 293 K and 2000 K, respectively).

The melting of the plate starts at . The time evolution of the accumulated melted volume of the solid material is illustrated in the graph of Figure 155.

Error creating thumbnail: File missing
Figure 155: Basic model. Time evolution of the melted volume of the steel plate.

From the graph one can deduce that the law is non-linear. This behavior is in part a consequence of the hypothesis of constant temperature for the corium. For this reason the melted particles of the solid elements take immediately the temperature of the corium that is greater than the melting temperature of the solid. A more reliable law would be obtained by considering variable the temperature of the corium and by using a higher heat capacity for the solid elements.

Another interesting point of the graph is represented by the increasing of the velocity of melting during the final instants of the simulation. The reason is that during this phase there is an increasing number of solid elements that become part of the fluid domain not due to the temperature criterion (see Figure 143) but due to the plasticity effect. In fact, in that moment of the analysis a huge part of the width of the plate has already melted and the thin remaining layer of the plate has to sustain all the weight of the corium and the melted steel. In this situation, that zone of the plate plastifies and it undergoes huge plastic deformations. A proof of this can be deduced from the graph of Figure 156.

Error creating thumbnail: File missing
Figure 156: Basic model. Time evolution of the temperature at the center of the lower side of the plate.

The evolution of the temperature in the plate bottom does not present that pronounced peak in the final part. This explains that the final increasing in the number of melted elements is not related directly with the effect of the temperature.

The collapsing of the plate occurs at . The hole created by the corium enlarges quickly due to the large deformations accumulated in the center of the plate.

This effect can be clearly seen in Figures 157 and 158, where several snapshots of this phase of the simulation are given.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 157: Basic model. Snapshots of the piercing of the melted plate by the corium (I/II).
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 158: Basic model. Snapshots of the piercing of the melted plate by the corium (II/II).

In Figure 159 the failure of the plate is represented in its deformed configuration enlarged by a factor of . From this representation it can be appreciated better the localization of the strains in the central part of the plate.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 159: Basic model. Deformed mesh of the plate during the last phase of the simulation (enlargement factor=).
In Figure 160 the plate deformation at the time is illustrated in a 3D-view with the same enlargement factor ().
Error creating thumbnail: File missing
Figure 160: Basic model. 3D-view of the deformed mesh of the plate at t=93.0 s (enlargement factor=).

In Figure 161 the XZ-view of the steel plate at the end of the analysis is given.

Error creating thumbnail: File missing
Figure 161: Basic model. Final configuration of the steel plate.

6.4 Detailed model

6.4.1 Problem data

The problem consists of a volume of corium at that leans against a steel shell and surrounds a rod that is located in the middle of the steel structure. In Figures 162 and 163 the initial geometry of the problem is illustrated. All the solid components of the problem have an initial temperature of . The objective of the analysis is to simulate the melting of the solid structure due to the heating caused by the corium.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
XZ view of the steel shell and corium
Figure 162: Detailed model. Initial geometry (dimensions in ).
Error creating thumbnail: File missing
Error creating thumbnail: File missing
XY view
Figure 163: Detailed model. 3D views of the initial geometry.

The shell is constrained in the same manner as the plate of the basic model Figure 152 while the rod is constrained at its top.

The material properties of the rod components, as well their dependency on the temperature, are the ones given in the specification and they are collected in Tables 15 and 16, respectively.


Table. 15 Material data for the rod components at T=293 K.
Inco. 600S 8430 216 0.33 14.86 444 0.282 1693
SUS310S 8030 200 0.33 13.31 502 0.242 1693

Material data for the rod components at T=293 K. MatDataRod


Table. 16 Material data for the rod components at different temperatures.
Inco. 600 216.7 215.7 209.9 198.1 197.1 187.3
SUS310S 200.1 - - - 170.6 157.9
Inco. 600 282 - - - - 220
SUS310S 242 - - - - -
Inco. 600 14.86 15.74 17.46 19.18 20.93 22.78
SUS310S 13.31 14.24 16.24 - 21.02 18.84
Inco. 600 444 465 486 507 528 553
SUS310S 502 - - - - -
Inco. 600 - 1.33E-5 1.39E-5 1.42E-5 1.46E-5 1.51E-5
SUS310S - 1.59E-5 - 1.62E-5 - 1.69E-5
Inco. 600 179.5 168.7 158.9 - -
SUS310S 145.1 131.4 114.7 - -
Inco. 600 216.7 179.5 75.5 41.2 -
SUS310S 146.1 133.4 84.3 44.1 21.6
Inco. 600 24.83 26.84 28.85 31.02 -
SUS310S 26.38 - 30.86 - -
Inco. 600 586 607 624
SUS310S - - - - -
Inco. 600 1.55E-5 1.60E-5 1.64E-5 1.67 E-5
SUS310S 1.73E-5 - - - 1.92E-5

Material data for the rod components at different temperatures. MatDataRodT

Concerning the plate and the corium, the material properties are the ones already introduced for the basic model (Figures 12, 13 and 14). As for the basic model, it has been assumed that the corium keeps the same temperature () for all the duration of the analysis.

6.4.2 Preliminary study

In the detailed model the steel structure is not plane as for the basic model, so the corium can not spread over the structure and fall by its edges. For this reason, in this case, it was not necessary to use the non-linear viscosity of Eq.(273).

However, with the proposed value of the viscosity , the corium should have a water-type behavior. In fact, the ratio of the corium is very similar to that of water:

(274)
(275)

It is well known that the Reynolds number depends on this ratio. Hence, the dynamics of the problem involving corium should be similar to that of the same problem involving water. With the proposed value of viscosity, the corium initially splashes over the shell, then starts oscillating until a hydrostatic state is reached. This is in contradiction with the description of corium as a 'highly viscous fluid' given in the specifications. Furthermore, from the computational point of view, for modeling this kind of problem a small time step increment should be used and the discretization of the upper surface of the plate should be fine. The union of these two requirements would increase highly the computational cost of the analysis.

For all these reasons, a higher value of viscosity has been considered. More specifically, has been used.

6.4.3 Numerical results

The finite element mesh used for the simulation consists of 396984 four-noded tetrahedra (71992 nodes). Both the fluid and the solid have been discretized using a non-structured mesh. In the central parts of the shell and the rod a finer mesh has been adopted. In particular, in those zones of the solid and in all the fluid domain the mesh has an average size of .

As for the basic model, the time step increment has been set depending on the phase of the simulation. In the first phase, when the corium is spreading on the shell and around the rod, the time step increment was used. In Figure 164, some snapshots of this phase are given. After the viscous fluid is almost at rest.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 164: Detailed model. Snapshots of the initial spreading of corium.

For the rest of the simulation the time step increment was increased to .

In Figure 165 some representative snapshots of the rest of the analysis are given.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 165: Detailed model. Snapshots of the final phase of the simulation. A hole is created in the rod and the corium passes through it.

The results show that the failure of the structure occurs in the rod. This can be clearly visualized in Figure 166 where only the structure is represented: the central part of the rod is completely melted, while no parts of the shell have reached the melting temperature of the material.

Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Error creating thumbnail: File missing
Figure 166: Detailed model. Snapshots of rod melting.

In the graph of Figure 167 the time evolution of the temperature on the external surfaces of the rod and the shell is plotted. This graph confirms that the rod is heating faster that the shell.

Error creating thumbnail: File missing
Figure 167: Detailed model. Time evolution of the temperature on the rod and shell external surfaces.

The melting of the rod starts at . The time evolution of the accumulated melted volume of the solid part of the domain is illustrated in the graph of Figure 155.

Error creating thumbnail: File missing
Figure 168: Detailed model. Time evolution of the melted volume.

Note that in this case the peak of melting velocity that occurs in the basic model (see Figure 155) does not happen. That is because in this problem the phase change of the solid is essentially governed by the temperature. The plastic zones are local and they do not create any plastic mechanisms, as it occurs for the basic problem (see Figure 159).

In Figure 169 the deformation of the solid structure at is illustrated in a 3D-view with an enlargement factor of .
Error creating thumbnail: File missing
Figure 169: Detailed model. 3D-view of the deformed configuration of the shell and the rod at t=15.0 s (enlargement factor=).

6.5 Summary and conclusions

In this chapter the potentialities of the PFEM have been tested for the solution of two very complex problems related to the melting of two structures in two specific situations of NPP severe accident. The simulations involved many critical phenomena for a finite element analysis, such as free-surface flows, thermal fluid-structure interaction, and phase change.

The objectives of the study have been reached for both the basic and detailed models proposed. In particular, the phenomena required by the contractor have been modeled successfully.

The analyses were run by using a slightly different model than the proposed one. Preliminary tests indicated that with the geometry and viscosity of the corium proposed initially by the contractor, the fluid would splash and fall down without melting the adjacent structure. For this reason, with the agreement of the contractor, some input data were modified. In particular, the corium was modeled using a higher value of the viscosity and, for the basic model, a smaller initial height was considered.

Due to the complexity of the analysis, some approximations and simplifications were accepted. The purpose of most of these was to reduce the computational cost of the analyses.

Among these assumptions, the most invasive ones were the use of the constant temperature for the corium and the simplified modeling of the different material properties of the solid melted volume. In particular, the constant temperature for the corium may have accelerated the melting of the structure. A more accurate simulation would be obtained by considering the corium to have a variable temperature and a higher heat capacity.

7 Conclusions and future lines of research

Chapter 7. Conclusions and future lines of research

The aim of this chapter is to summarize all the work done for this thesis highlighting its innovative points and main contributions. The thesis opens new perspectives to several developments. The future lines of work are given in the last section of this chapter.

7.1 Contributions

The objective of this thesis was the derivation and the implementation in a code of a unified formulation for fluid and solid mechanics, fluid-structure interaction and thermal coupled problems using the PFEM.

The finite element procedure was derived starting from the Velocity formulation presented in Chapter 2. In the same chapter, the mixed Velocity formulation was obtained exploiting the linearized form of the Velocity formulation. The mixed Velocity-Pressure scheme is based on a two-step Gauss-Seidel procedure where first the linear momentum equations are solved for the velocity increments and then the continuity equation is solved for the pressure in the updated configuration. Linear interpolation was used for both the velocity and the pressure fields. Both velocity-based finite element Lagrangian procedures were derived first for a general compressible material and then particularized for the hypoelastic model. The hypoelastic solid elements generated from the Velocity formulation and the mixed Velocity-Pressure formulations were called V and VP elements, respectively. The V and the VP elements were validated for several benchmark problems for elastoplastic compressible structures. It has been shown that both elements are convergent for all the numerical examples analyzed.

In Chapter 3 the Unified Stabilized formulation for quasi-incompressible materials was derived. This numerical procedure essentially consists of the mixed Velocity-Pressure formulation derived in Chapter 3 where for the continuity equation the FIC stabilized form is used. The FIC stabilization was derived for the case of quasi-incompressible fluids and was extended also to hypoelastic quasi-incompressible solids. The stabilized element for quasi-incompressible solids was called VPS element. The solution schemes for solving quasi-incompressible Newtonian fluids and hypoelastic solids with the stabilized mixed Velocity-Pressure formulation were given and explained in detail. The entire Section 3.4 was devoted to the analysis of free surface fluids with the Unified stabilized formulation. First the Particle Finite Element Method (PFEM) was explained analyzing its advantages and disadvantages. An innovative strategy for modeling slip conditions in a Lagrangian way has also described. The excellent mass preservation properties of the PFEM-FIC stabilized formulation by solving a variety of 2D and 3D free surface flow problems involving surface waves, water splashing, violent impact of flows with containment walls and mixing of fluids. Also the conditioning of the scheme was studied and the effect of the bulk modulus on the numerical scheme highlighted. It has been shown that using a scaled value for the bulk modulus, pseudo bulk modulus, in the linear momentum tangent matrix improves the conditioning and the global convergence of the linear system. A simple and efficient strategy for calibrating the optimum value for the pseudo bulk modulus was also given. The strategy was successfully validated for two benchmark problems for free surface flow. In the last section of Chapter 3 the unified stabilized formulation for fluids and solids at the incompressible limit was validated comparing the numerical results to both experimental tests and numerical results from other formulations. It was shown that the method is convergent to the expected solution. Specific attention was given to the study of the boundary conditions highlighting the beneficious effect of using slip conditions for inviscid fluids and coarse meshes.

Chapter 4 was devoted to the application of the Unified Stabilized to FSI problems. It was shown that this operation requires a small implementation effort. It is only required to assemble properly the global linear system and to detect the interface exploiting the capability of the PFEM for detecting the contours. It was shown that, depending on the specific need, one may choose any of the V, the VP and the VPS elements for modeling hypoelastic solids. The numerical solution given by the proposed scheme was validated with analytical solutions and numerical results of other formulations for free surface FSI problems.

In Chapter 5 the Unified formulation was used for solving coupled thermal-mechanical problem. The coupling was ensured via a staggered scheme and for the temperature field the same linear shape functions of the velocity and pressure fields were used. The numerical method was tested with several solid and fluid dynamics problems involving the temperature. The algorithm for the phase change modeling was explained and an explicative numerical example was given.

In Chapter 6 the Unified stabilized thermal coupled formulation has been tested with an industrial problem. The analysis concerned the melting of two structures in two specific situations of Nuclear Power Plant (NPP) accident. The simulations involved many critical issues for a finite element analysis, such as free-surface flows, thermal fluid-structure interaction and phase change. Because of the complexity of the analysis, some approximations and simplifications were accepted by the contractor. The project ended successfully and all the phenomena required by the contractor in the specifications were modeled.

The innovative points and the contributions of this thesis can be summarized by the following list:

7.2 Lines for future work

This section outlines the possible lines of research opened by this thesis. The following enumeration summarizes most of these

toc

toc

dummy ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––

Bibliography

BIBLIOGRAPHY

Back to Top

Document information

Published on 01/01/2011

Licence: CC BY-NC-SA license

Document Score

0

Views 0
Recommendations 0

Share this document

claim authorship

Are you one of the authors of this document?